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CHAPTER  1 
INTRODUCTION 


There  has  been  increased  interest  in  shear  wave  exploration  in  the  past 
decade.  A  thorough  knowledge  of  the  weathered  zone  shear  wave  velocities 
can  be  used  for  static  corrections  in  surface  seismic  applications  (Mari,  1984).  If 
the  information  on  velocity  structures  extends  deep  enough,  then  the  determined 
shear  structure  can  be  used  directly  in  exploration.  Shear  and  compressional 
velocities  can  be  related  to  elastic  parameters  (Waters,  1981): 


(1.1) 


|3  =  /ii  s=  shear  velocity, 

Ve 


K  +  4  / 


^  =  compressional  velocity. 


where  K  is  the  bulk  modulus  of  elasticity,  q  is  the  bulk  density,  and  p  is  the 
modulus  of  rigidity.  Another  common  way  of  looking  at  material  properties  is 
through  the  use  of  Poisson's  ratio  which  can  be  expressed  in  terms  of  shear  and 
compressional  velocities: 

p  p 

, .  transverse  strain  -  23*-  „  .  ,  *• 

(1.2)  a- - - - - — =  Poisson  s  ratio. 

longitudinal  strain  2a^-2|3^ 
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As  a  result,  a  knowledge  of  S-wave  velocities  can  be  used  to  determine  elastic 
moduli  of  earth  materials  at  shallow  depths.  This  information  has  substantial 
engineering  application  (Kitsunczaki,  1980).  For  earthquake  engineering,  S- 
wave  velocities  are  used  for  evaluation  of  transmission  characteristics  of  the 
ground.  This  information  in  turn  leads  to  knowledge  of  ground  vibrations  during 
earthquakes,  which  can  be  used  to  minimize  earthquake  damage.  Elastic 
moduli  determination  through  shear  velocities  is  also  used  for  site  evaluation  in 
civil  engineering. 

Modeling  of  shallow  layers  or  the  weathered  zone  has  historically  been 
accomplished  with  refraction  surveys.  Miller  and  Steeples  (1990)  have  utilized 
reflection  techniques  to  yield  structure  and  compressional  velocity  as  shallow  as 
30  meters.  The  refraction  method  works  quite  well  for  determining  compres¬ 
sional  wave  velocities  and  layer  thicknesses.  Defining  the  S-wave  velocity  struc¬ 
ture  is  not  generally  as  easily  accomplished.  Since  P-wave  velocities  are  always 
faster  than  shear  wave  velocities,  S-wave  arrivals  are  always  contaminated  by 
compressional  arrivals.  An  S-wave  refraction  technique  could  only  work  if  a 
source  generated  shear  waves  without  compressional  waves  or  if  the  two  wave- 
fields  could  be  easily  separated.  It  has  historically  been  very  difficult  to  generate 
shear  waves  without  creating  an  interfering  compressional  disturbance.  Picking 
the  first  arrival  time  on  a  shear  wave  thus  very  difficult  because  of  this  P-wave  in¬ 
terference.  Additional  interference  is  caused  by  large  amplitude  surface  waves 
or  ground  roll,  which  has  a  velocity  similar  to  the  shear  velocity  of  the  surface 
layer.  This  study  will  include  data  acquired  with  a  new,  clean  SH  energy  source. 
Refraction  processing  of  both  P  and  SH  data  will  be  employed  prior  to  an  inter¬ 
pretation  of  the  surface  wave  arrivals. 
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There  is  imaging  information  applicable  to  the  near  surface  structure 
throughout  the  wavetrain  of  short  array  seismic  data  sets.  The  body  wave  por¬ 
tion  includes  reflection  and  refraction  information.  Most  body  wave  techniques 
require  the  minimization  of  the  surface  wave  portion  of  the  wavetrain.  The  sur¬ 
face  waves  can  be  used  separately  to  determine  near  surface  material  models. 
In  this  study,  the  body  wave  models  will  be  compared  to  surface  wave  models. 
One  technique  that  has  been  used  successfully  to  determine  shear  velocity 
structures  is  the  interpretation  of  surface  wave  dispersion  curves.  The  technique 
involves  measuring  the  surface  wave  velocity  as  a  function  of  frequency  or  peri¬ 
od.  The  measured  surface  wave  dispersion  is  then  related  to  shear  velocity 
structure.  A  shear  wave  source  is  used  to  generate  Love  waves  and  a  com- 
pressional  source  is  used  to  generate  Rayleigh  waves.  Both  phase  and  group 
velocity  dispersion  curves  from  Love  and  Rayleigh  waves  will  be  analyzed.  Re¬ 
lating  dispersion  curves  to  velocity  structure  is  done  with  an  iterative  forward 
modeling  technique  or  one  of  a  number  of  inversion  techniques.  Both  ap¬ 
proaches  will  be  utilized  in  this  study. 


1.  Surface  Waves 


Since  the  new  aspect  of  this  study  is  the  interpretation  of  the  surface 
wave  portion  of  the  wavefield,  some  fundamental  concepts  that  deal  with  surface 
waves  and  surface  wave  propagation  are  introduced.  Surface  wave  effects  are 
confined  closely  to  the  free  surface  boundary.  Rayleigh  waves  are  a  type  of  sur¬ 
face  wave  that  are  trapped  near  a  free  plane  boundary.  Rayleigh  waves  consist 
of  a  mixture  of  shear  (SV)  and  compressional  energy.  For  the  theoretical  case 
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of  a  half  space,  the  velocity  along  the  surface  is  independent  of  frequency  and 
depends  on  the  shear  and  compressional  velocities  or  equivalently,  Poisson's 
ratio.  For  a  Poisson's  ratio  of  1/4,  the  speed  of  the  Rayleigh  wave  is  .9194  times 
the  shear  velocity.  The  particle  motion  at  the  surface  is  retrograde  elliptical.  The 
amplitude  of  motion  decays  away  from  the  free  surface  as: 

(1 .3)  Amp  =  exp(-r  Z  2:i/X), 

where  r  is  a  factor  related  to  velocity  and  Z  is  depth.  The  amplitude  decays 
more  rapidly  for  small  wavelengths  or  high  frequencies.  As  a  result,  long  wave¬ 
lengths  penetrate  deeper  below  the  boundary.  In  most  normal  weathered  situ¬ 
ations,  velocities  increase  with  depth.  This  velocity  profile  leads  to  the  situation 
of  normal  dispersion  in  phase  velocity,  where  long  wavelengths  have  a  greater 
phase  velocity  than  short  wavelengths. 

Love  waves  are  large  amplitude  trapped  surface  waves  whose  presence 
depends  on  near  surface  layering.  Their  presence  makes  the  interpretation  of 
shear  wave  (SH)  reflection  records  difficult.  A  simple  model  for  a  Love  wave  can 
be  illustrated  by  considering  an  SH  wave  propagating  from  a  source,  S,  in  a 
layer  over  a  half  space  (Figure  1).  The  wave  front  velocity  or  phase  propagation 
Cs,  along  the  surface  can  be  shown  geometrically  to  be: 

(1 .4)  Cs  =  (3i/cos(n/2-e), 

where  (3i  is  the  shear  velocity  in  layer  1.  TTie  critical  angle  Sq  is: 
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(1.5)  0Q  =  arcsin  Pi/P2, 

where  (32  is  the  shear  velocity  in  the  halfspace.  The  larger  the  contrast  for  |3,  the 
smaller  Sc  becomes.  For  any  wave  propagating  with  an  angle  equal  to  or 
greater  than  the  critical  angle,  the  wave  will  be  trapped  between  the  layer  and 
the  surface  interfaces.  No  mode  conversions  will  occur  because  only  SH  waves 
are  considered.  A  phase  shift  dependent  on  the  incident  angle  will  occur  at  the 
layer  interface,  but  there  will  be  no  phase  change  for  SH  at  the  free  surface. 
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Figure  1 .  Geometric  analysis  of  Love  wave  generation  in  layer 
over  a  half  space  (After  Waters,  1981). 


For  waves  of  a  constant  frequency,  there  are  points  on  the  same  wave- 
front,  such  as  E  and  F  (Figure  1),  that  constructively  interfere  if  the  distance 
EBCF  is  an  integral  multiple  of  the  wavelength.  If  a  constant  phase  change  at 
the  layer  interface  is  assumed,  it  can  be  argued  geometrically  that  the  construc¬ 
tive  modes  will  be  dependent  only  on  the  wavelength  and  the  angle  of  incidence. 
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As  the  wavelength  or  frequency  changes,  then  6  must  change  to  maintain  con¬ 
structive  interference.  Since  equation  (1 .4)  has  already  demonstrated  that  the 
phase  velocity  is  dependent  on  0,  it  follows  that  the  phase  velocity  of  this  con¬ 
structive  interference  changes  with  frequency.  This  frequency  dependence  of 
velocity,  known  as  dispersion,  causes  the  wave  packet  to  spread  out  as  it  propa¬ 
gates.  This  dispersion  is  the  reason  that  phase  velocity  and  group  velocity  are 
not  identical. 


2.  Phase  and  Group  Velocity 

The  distinction  between  phase  and  group  velocity  is  fundamental  to  the 
dispersion  techniques  utilized  in  this  thesis.  Velocity  is  simply  the  rate  of  propa¬ 
gation  in  a  medium.  Velocity,  V,  is  proportional  to  the  distance  traveled  D,  di¬ 
vided  by  the  time  of  travel  At: 

(1.6)  y  =  —. 

At 

Phase  velocity  is  associated  with  the  propagation  velocity  of  a  particular  feature 
of  the  waveform,  such  as  a  peak  or  a  trough  (Waters,  1981).  The  term  phase 
velocity  refers  to  the  effective  velocity  along  the  ground  surface  of  a  particular 
frequency.  If  observations  are  made  at  points  sufficiently  close  together,  the 
phase  velocity  C  is  simply  the  distance,  AD,  between  observation  points  divided 
by  the  change  in  time  of  arrival  of  a  particular  phase.  At: 


(1.7) 


C  =  AD/At  . 
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Phase  velocity  can  also  be  thought  of  as  the  change  in  phase  angle  0,  meas¬ 
ured  tor  the  frequency  t,  over  time,  which  would  be  proportional  to  angular  fre¬ 
quency  oj.  This  linear  phase  shift  is  proportional  to  a  time  delay: 

(1.8)  At  =  0{o))/a). 

This  results  in: 

(1.9)  C  =  D/At  =  ajD/0  =  2:ifD/0. 

It  has  been  previously  shown  that  for  some  surface  waves,  the  phase  velocity  is 
dependent  upon  frequency. 

Group  velocity  can  be  defined  as  the  velocity  with  which  the  energy  asso¬ 
ciated  with  a  narrow  band  of  frequencies  has  traveled.  Group  velocity  is  a  func¬ 
tion  of  frequency,  the  velocity  of  the  media,  and  the  thickness  of  the  layering. 
The  group  velocity  U  for  a  given  period  T,  is  related  to  the  phase  velocity,  C: 

(1.10)  U(T)  =  C  +  k^, 

dk 

where  k  «  wave  number  =  2:t/CT. 

In  rts  simplest  form,  group  velocity  is  inversely  proportional  to  the  time  for 
a  delta  function  signal  to  travel  from  a  source  to  a  receiver.  If  there  is  no  disper¬ 
sion  or  variation  in  velocity  with  frequency  or  period,  then  the  group  and  phase 
velocities  are  identical.  Even  with  dispersion,  the  group  velocity  can  be  deter- 
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mined  with  only  one  seismometer,  since  the  group  velocities  can  be  calculated 
by  the  simple  equation; 


(1.11)  U(T)  =  ^. 

At 

Figure  2  displays  group  and  phase  velocities  as  a  function  of  frequency  for 
Rayleigh  waves  generated  from  a  plane  layered  model.  The  model  was  derived 
from  P  and  SH  wave  refraction  analysis  of  data  acquired  for  this  study.  The  Q 
values  utilized  were  those  from  the  Grant  (1988)  model  (Table  1).  The  top  two 
curves  are  the  first  higher  mode.  The  higher  velocity  curve  in  each  pair  is  the 
phase  velocity,  while  the  lower  curve  Is  the  group  velocity.  Since  the  group  ve¬ 
locity  is  related  by  a  derivative  to  the  phase  velocity  (equation  1.10),  the  resulting 
group  velocity  curve  is  less  smooth  than  the  corresponding  phase  velocity  curve. 

3.  Dispersion  Analysis  Using  Paper  Records 
One  of  the  earliest  techniques  for  determining  group  velocity  of  dispersed  sur¬ 
face  waves  was  known  as  the  peak  and  trough  technique.  Press  (1956)  docu¬ 
ments  how  this  procedure  was  employed.  He  would  make  use  of  the  paper  re¬ 
cords  of  the  Rayleigh  wave  primary  mode  arrival  from  one  earthquake  on  a  sin¬ 
gle  seismic  recording  station.  The  individual  cycles  were  measured  to  yield  the 
period  T,  to  be  used  with  the  known  distance  D  and  measured  At  or  travel  time 
from  the  epicenter.  This  process  could  be  repeated  for  a  series  of  periods  to  ob¬ 
tain  the  group  velocity  dispersion  curve.  This  group  velocity  would  be  a  function 


VRLOCITY  (M/SE£C) 


FREQUENCY  (HZ) 


Figure  2.  Rayleigh  wave  dispersion  curves  from  a  plane  layered 
synthetic  model.  Curves  from  top  to  bottom  are:  RPH1 ,  Rayleigh 
wave  phase  velocity,  higher  mode  1.  RGH1,  Rayleigh  wave  group 
velocity,  higher  mode  1.  RPF,  Rayleigh  wave  phase  velocity, 
fundamental  mode.  RGF,  Rayleigh  wave  group  velocity, 
fundamental  mode. 
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of  the  average  velocities  of  the  crust  between  the  point  of  origin  for  the  earth¬ 
quake  and  the  recording  station.  To  determine  local  variations,  it  is  necessary  to 
use  the  phase  velocity.  At  least  two  receivers  must  be  utilized  to  calculate  phase 
velocity.  The  peak  and  trough  technique  for  phase  velocity  determination  was 
also  documented  by  Press  (1956).  TTiis  method  can  best  be  understood  by  us¬ 
ing  two  recordings  along  a  great  circle  path  from  the  epicenter  of  an  earthquake. 
If  the  two  stations  are  within  one  wavelength  of  each  other,  the  two  seismograms 
of  the  Rayleigh  wave  portion  can  be  aligned  in  time  and  then  carefully  compared. 
The  peaks  and  troughs  that  are  the  same  can  be  compared  for  time  differences. 
The  time  difference  is  At  and  the  distance  D  can  be  measured  along  the  great 
circle  path.  The  period  T  is  measured  as  before,  by  the  length  of  the  cycle.  The 
actual  implementation,  in  practice,  is  slightly  more  complicated  since  the  earth¬ 
quakes  and  seismometers  do  not  normally  fall  conveniently  on  a  great  circle 
path.  Three  stations  and  some  simple  geometry  are  used  to  make  the  actual 
phase  velocity  calculations.  The  final  dispersion  curve  can  be  used  as  before  to 
make  inferences  about  the  crustal  structure,  but  the  horizontal  extent  is  much 
smaller  than  before. 

The  analyi  is  of  the  dispersion  results  was  done  by  computing  theoretical 
phase  velocity  curves  from  measured  group  velocities  using  equation  (1.10). 
These  group  velocities  were  interpreted  to  be  the  results  of  the  average  structure 
and  average  elastic  constants  over  a  long  distance.  The  measured  phase  ve¬ 
locities  were  then  compared  to  these  predicted  curves.  Any  differences  were  in¬ 
terpreted  to  be  due  to  local  variations  in  crustal  thickness. 
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4.  Dispersion  Analysis  Using  Moving  Window  Technique 

The  previous  technique  was  utilized  before  the  widespread  availability  of 
digital  computers  and  the  advent  of  the  FFT.  The  peak  and  trough  technique 
would  not  work  well  under  conditions  of  low  signal  to  noise  ratio  or  if  the  signal 
was  otherwise  obscured.  Microseisms  or  long  period  drifts  may  obscure  the  sig¬ 
nal.  Additionally,  several  dispersive  modes  or  several  branches  of  the  same 
mode  may  appear  simultaneously  at  the  recording  site. 

The  moving  window  analysis  (Landisman  et  al.,  1969)  can  be  used  to 
measure  group  velocities  in  recordings  of  multi-mode  propagation.  The  time  and 
distance  measured  from  the  earthquake  epicenter  yield  a  group  velocity  for  a 
point  on  the  time  series.  A  section  of  the  time  series  centered  around  this  point 
or  group  velocity  is  extracted  by  windowing  with  a  boxcar  in  the  time  domain. 
The  extracted  data  is  then  windowed  again  with  a  symmetrical  weighting  function 
such  as  cos2  in  the  time  domain.  Windowing  the  second  time  reduces  the 
ringing  or  Gibb's  phenomena  in  the  frequency  domain.  The  symmetrical  nature 
of  the  weighting  function  emphasizes  the  center  of  the  window  which  most 
closely  corresponds  to  the  current  time  or  group  velocity  of  interest.  The  next 
stage  of  moving  window  analysis  involves  the  Fourier  transformation  of  the  win¬ 
dowed  data  and  the  presentation  of  instantaneous  amplitude  in  dB  versus  peri¬ 
od.  This  process  is  continued  over  time  or  group  velocity.  The  window  is  of 
variable  length  in  order  to  maintain  the  same  frequency  resolution  at  all  periods. 
After  complete  analysis  a  plot  is  generated  on  which  equal  energy  contours  are 
drawn  as  shown  in  Figure  3.  Figure  3  contains  moving  window  analysis  of  fun¬ 
damental  mode  Rayleigh  wave  data  from  a  study  by  Simila  (1982)  in  the  same 
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Figure  3.  Moving  window  analysis  (Simila,  1982).  Group  velocity 
analysis  of  Rayleigh  wave  data  at  McCormick  Ranch.  T  is  period 
milliseconds,  U  is  group  velocity  in  m/s. 
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area  as  this  thesis.  The  hand  drawn  circles  through  the  maximum  energy  at 
each  period  and  velocity  define  the  group  velocity  or  dispersion  curve, 

5,  Inversion  Fundamentals 

Forward  modeling  is  defined  (Menke,  1 989)  as  the  process  of  predicting 
the  results  of  measurements  on  the  basis  of  some  general  principle  or  model 
and  a  set  of  specific  conditions  relevant  to  the  problem  at  hand.  Inverse  theory 
tackles  the  reverse  problem,  starting  with  measured  data  and  a  general  principle 
or  model,  it  determines  estimates  of  the  model  parameters.  An  example  of  in¬ 
verse  theory  or  inversion  would  be  determining  a  velocity  from  a  given  set  of 
time  offset  pairs.  For  a  presumed  constant  velocity,  the  problem  would  be  to  fit  a 
straight  line  through  a  series  of  points.  With  real  data,  in  general,  the  points 
would  not  all  fall  exactly  on  a  line.  Part  of  inverse  theory  would  involve  deriving  a 
best  fit  of  a  model  to  the  data. 

The  physical  properties  that  can  be  dealt  with  using  inversion  techniques 
fall  into  two  general  classes.  The  first  category  would  be  properties  that  can  be 
described  by  discrete  parameters  and  the  second  class  would  be  those  that  are 
described  by  continuous  functions.  This  paper  will  deal  with  discrete  inverse 
theory  or  parameters  that  are  truly  discrete  or  that  can  be  approximated  as  dis¬ 
crete.  The  parameters  that  are  determined  in  this  thesis  are  assumed  to  have  a 
finite  quantity  of  numerical  values.  The  earth  will  be  assumed  to  be  composed  of 
several  discrete  layers  tha^  can  be  represented  by  fundamental  parameters,  in¬ 
cluding  velocity,  density  and  Q  or  attenuation.  The  data  used  to  determine  these 
parameters  will  be  the  surface  wave  phase  and  group  velocity  dispersion  curves 
Introduced  earlier. 
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One  of  the  most  easily  understood  forms  of  inverse  problems  is  the  linear 
equation: 

(1.12)  Gm  =  d. 

where  d  is  a  vector  of  the  data,  m  is  a  vector  or  matrix  of  the  model  parameters 
being  estimated  and  G  is  the  data  kernel.  This  equation  is  at  the  heart  of  most 
discrete  inverse  theory.  Inversion  Is  generally  an  iterative  technique  where  this 
equation  can  be  worked  backwards  based  on  the  estimated  model  parameters 
most  to  yield  the  predicted  data  dP*'®  =  G  m®s^.  The  frequently  used  least 
squares  solution  is  a  technique  to  pick  model  parameters  that  force  the  predicted 
data  to  be  as  close  as  possible  to  the  observed  data.  An  error  or  misfit  for  each 
observation  can  be  defined  as: 

(1-13)  e.  =  -  dr*- 

The  best  fit  line  is  then  the  one  with  the  smallest  summed  squared  error  E,  de¬ 
fined  as: 


(1.14) 


This  technique  is  common  for  deriving  a  best  fit  line  to  a  series  of  points.  A  non¬ 
linear  problem  can  be  solved  by  approximating  it  as  a  linear  problem.  The 
model  can  be  slightly  perturbed,  and  the  comparison  of  the  predicted  data  to  that 


observed  can  be  checked.  An  estimate  of  the  improved  fit  can  be  made  as  well 
as  iterative  changes  to  the  model.  As  with  most  iterative  schemes,  large  pertur¬ 
bations  can  lead  to  either  overshoot  or  undershoot.  The  perturbation  of  the 
model  can  be  done  using  a  partial  derivative  approach,  which  is  effectively  the 
relationship  between  the  changes  of  two  parameters.  The  inversion  techniques 
used  in  this  thesis  are  based  upon  the  principles  documented  in  this  section. 

B.  Previous  Work 

The  basic  technique  of  using  surface  wave  dispersion  to  model  shear  ve¬ 
locity  structures  has  been  used  by  several  groups.  Mooney  and  Bolt  (1 966)  cal¬ 
culated  theoretical  group  and  phase  velocity  dispersion  curves  for  Rayleigh 
waves  under  varying  conditions  of  shear  velocity  contrast  for  a  layer  over  a  half 
space.  Included  in  their  numerical  exercise  was  a  model  for  a  shallow  alluvium 
site  that  included  the  calculation  of  the  fundamental  as  well  as  the  first  and  sec¬ 
ond  higher  modes.  The  results  indicated  that  periods  of  up  to  2  seconds  would 
be  required  to  image  1 00  meters  of  alluvium  or  a  0.5  second  period  to  image  25 
meters. 

Mari  (1984)  used  Love  wave  phase  velocity  measurements  to  obtain 
static  corrections.  He  generated  Love  waves  with  the  Marthor  source,  which  is  a 
high  energy  shear  impulsive  source.  TTiis  study  utilized  exploration  seismology 
phones,  as  well  as  normal  exploration  seismology  source  and  receiver  spreads. 
The  source-receiver  spacings  ranged  from  250  to  490  meters.  The  results  were 
generated  using  frequencies  from  5  to  20  Hz  to  yield  surface  layer  thicknesses  of 
6  to  10  meters  and  phase  velocities  of  140  to  210  m/s.  This  study  built  on  the 
previous  numerical  analysis  of  Mooney  and  Bolt  (1 966).  Mari  assumed  a  layer 
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over  a  half  space  model,  and  interpreted  his  slowest  velocity  as  the  shear  veloc¬ 
ity  tor  the  layer  and  the  fastest  velocity  as  the  velocity  for  the  half  space  or  sub¬ 
stratum. 

Cherry  et  al.  (1979)  used  Rayleigh  wave  group  velocities  to  create  a 
shear  model.  This  study  used  a  weight  drop  source  and  a  linear  array  witl=r  re- 
ceivers  spaced  500  meters  apart  over  the  interval  1  to  5.5  km.  Surface  wave 
velocities  ranged  from  400  to  700  m/s.  Frequencies  between  2  and  20  Hz  were 
used  to  invert  for  shear  structure  to  500  meters.  Mokhtar  et  al.  (1988)  used 
phase  and  group  velocity  dispersion  to  invert  for  shallow  crustal  structure. 
Rayleigh  wave  frequencies  between  1  and  20  Hz  were  used  to  image  the  top  1 
to  2  km  of  the  crust.  The  source  receiver  spacings  were  55  to  80  km  and  yielded 
phase  velocities  in  the  range  of  2.8  to  3.4  km/s.  Velocities  were  determined  with 
linear  inversion  software  developed  by  the  authors.  The  same  software  is  util¬ 
ized  in  this  thesis. 

Surface  wave  investigations  have  been  completed  over  a  wide  range  of 
length  and  frequency  scales.  Press  (1956)  studied  crustal  structure  with  thick¬ 
nesses  from  15  to  50  km  using  Rayleigh  wave  phase  velocity  dispersion  from 
periods  of  10  to  40  seconds  (wavelengths  from  50  to  100  km).  The  Press  (1956) 
study  can  be  contrasted  to  the  study  of  Stokoe  and  Nazarian  (1984),  \who  also 
used  Rayleigh  wave  phase  velocity  dispersion.  This  work  was  to  determine 
shear  velocities  in  concrete  runways  approximately  ten  inches  thick.  This  re¬ 
markable  range  of  scale  can  be  achieved  by  utilizing  data  of  different  wave¬ 
lengths. 


This  study  will  utilize  seismic  data  acquired  in  a  manner  similar  to  a  nor¬ 
mal  refraction  survey.  The  test  site  is  the  McCormick  Ranch  area,  which  is  outh 
of  Albuquerque,  New  Mexico.  This  area  is  a  relatively  flat  region  at  the  base  of 
the  Sandia  Mountains.  The  subsurface  geology  consists  of  dry  alluvium  down  to 
the  water  table  at  approximately  75  meters.  Caliche  is  present  in  the  test  site 
and  has  been  exposed  during  some  high  explosive  test  shots.  The  layers  of  al¬ 
luvium  are  assumed  to  be  horizontally  stratified  and  laterally  homogenous. 
Some  of  the  previous  work  on  this  test  site  includes  a  study  by  Grant  (1988), 
who  performed  a  source  characterization  study  in  the  area  using  small  explo¬ 
sives  sources  which  required  some  ^nowiedge  of  the  near  surface  structure. 
Additional  studies  performed  in  the  area  which  characterized  the  shallow  struc¬ 
ture  includes  Stump  and  Reinke  (1982)  and  Reinke  and  Stump  (1988). 
Bogaards  (1989)  used  model  based  spectral  analysis  techniques  to  generate 
dispersion  curves  for  surface  waves  utilizing  similar  acquisition  techniques  to  the 
current  study.  Two  previous  models  for  the  McCormick  ranch  study  area  are 
summarized  in  Table  1. 

Simita  (1982)  used  Rayleigh  wave  group  velocities  to  model  velocity  struc¬ 
tures  in  shallow  (less  than  25  meters)  stmctures  at  McCormick  Ranch.  He  util¬ 
ized  explosive  sources  to  generate  the  Rayleigh  waves.  The  dispersion  analysis 
incorporated  single  phone  processing  at  distances  ranging  from  6.55  to  228.6 
meters.  The  processing  included  the  moving  window  analysis  technique  for  fre¬ 
quencies  from  5  to  25  Hz  for  the  fundamental  mode  and  25  to  50  Hz  for  the  first 
higher  mode.  This  thesis  will  extend  Simiia's  work.  New  processing  schemes 
are  available  to  analyze  arrays  of  data  for  both  group  and  phase  velocity.  New 
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equipment  is  now  available  that  makes  it  more  practical  and  economical  to  ac¬ 
quire  an  array  of  data,  as  opposed  to  the  single  phone  acquisition  used  in  the 
Simila  study.  A  new  shear  source  makes  it  possible  to  generate  shear  data  and 
a  resulting  Love  wave.  The  final  extension  of  this  work,  as  proposed  by  Simila, 
will  be  the  use  of  inversion  techniques,  instead  of  relying  on  a  trial  and  error 
forward  modeling  method. 

D.  Overview  of  the  Thesis 

The  seismic  data  acquisition  for  this  study  will  be  discussed  at  length  in 
Chapter  2.  The  processing  techniques  for  deriving  dispersion  curves  and  analy¬ 
sis  of  the  results  are  contained  in  Chapter  3.  Chapter  4  covers  inversion  of  the 
dispersion  data  and  retraction  data  for  velocity  structures.  Chapter  5  summa¬ 
rizes  this  work  and  includes  recommendations  for  further  study. 


Table  1 

PREVIOUS  MCCORMICK  RANCH  SITE  MODELS 


Grant  (1988)  Model 


a(m/s) 

p(m/s) 

(gm/cm3 

Qa 

Depth  (m) 

270 

120 

1.8 

10 

4 

0.0 

670 

230 

1.9 

50 

22 

3.0 

Half 

space 

930 

360 

2.0 

100 

44 

19.0 

Stump-Reinke  1982  M 

odel 

a(m/s) 

p(m/s) 

jgm/cm3 

Qa 

Depth  (m) 

lESDH 

366 

244 

1.8 

0.0 

671 

366 

1.9 

3.4 

lESISBi 

823 

366 

2.0 

13.4 

Half 

space 

1128 

610 

2-1 

24.0 

Where:  a  =  P-wave  velocity  in  m/sec 


p  =  S-wave  velocity  in  m/sec 

Q  =  density  in  gm/cm^ 
s  P-wave  attenuation 

Qj3  =  S-wave  attenuation 


CHAPTER  2 

FIELD  DATA  ACQUISITION 


The  seismic  data  utilized  in  this  study  was  acquired  by  Pearson  and 
Stump  of  SMU.  The  acquisition  hardware  was  fully  documented  by  Stump  et  al. 
(1991).  The  intent  for  this  study  was  to  generate  two  types  of  surface  waves, 
both  Rayleigh  and  Love  waves.  The  seismic  data  were  acquired  in  two  parts 
with  a  different  source  in  each  case.  The  P-SV  data  set  was  acquired  utilizing  a 
commercially  available  source,  the  Betsy  Seisgun*  (Figure  4).  The  Betsy  is  ef¬ 
fectively  a  shotgun  with  the  barrel  pointing  down  at  the  ground.  The  source  is  an 
eight  gauge  shotgun  shell  with  a  single  solid  lead  or  iron  slug.  A  tire  at  the  base 
of  the  source  is  used  to  muffle  airblast.  There  are  a  number  of  advantages  to 
using  the  Betsy  for  this  type  of  study.  The  source  is  relatively  portable,  simple 
and  inexpensive  to  operate,  and  is  capable  of  generating  a  repeatable  signature. 
This  source  is  relatively  low  energy,  but  is  capable  of  generating  sufficient  en¬ 
ergy  for  the  100  meters  of  maximum  source-receiver  offset  acquired  for  this 
study.  Frequency  content  of  this  source  is  in  excess  of  200  Hz,  although  the 
data  acquired  at  McCormick  Ranch  had  somewhat  less  bandwidth  due  to  the 
highly  attenuative  nature  of  the  media.  The  Betsy  generates  a  good  compres- 
sional  wavefield,  and  for  this  study  a  significant  Rayleigh  wave. 


*Note:  Betsy  Seisgun  is  a  registered  trademark  of  Betsy  Seisgun,  Inc. 
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The  second  source  was  the  SMU  fabricated  SWIG,  which  stands  for 
Shear  Wave  Impulsive  Generator.  The  SWIG  is  a  modification  of  a  source  de¬ 
scribed  by  Liu  et  al.  (1988).  The  SWIG  (Figure  5)  is  a  pneumatic  ram  mounted 
between  two  channels  that  is  weighed  down  by  a  pickup  truck.  The  SWIG  is  ori¬ 
ented  perpendicular  to  the  source-receiver  plane.  One  side  of  the  piston  is 
evacuated  to  the  outside  air  and  the  high  pressure  side  drives  the  mass  horizon¬ 
tally  until  it  impacts  an  anvil  at  the  end  of  the  SWIG.  This  impact  produces  a 
wavefield  with  particle  motion  perpendicular  to  the  direction  of  propagation  from 
the  source  to  the  geophones. 

The  SWIG  provides  a  very  clean  SH  wavefield.  In  the  test  site  for  this 
study,  the  source  generated  a  very  distinct  Love  wave  along  with  the  shear 
wavefield.  Figures  6  and  7  show  the  SH  and  P  first  arrivals  respectively.  It  can 
be  seen  that  there  is  no  arrival  corresponding  to  the  P  first  arrival  from  the  Betsy 
source  occurring  on  the  SWIG  data  set.  This  uncontaminated  shear  wavefield 
will  be  used  in  refraction  calculations  to  determine  shear  velocity  structure  as 
well  as  for  Love  wave  dispersion  analysis. 

The  data  was  acquired  sixty  channels  at  a  time.  The  recording  instruments  were 
six  channel  Refraction  Technology  data  acquisition  subsystems  (DAS),  time 
synchronized  through  a  radio  telemetty  system.  In  this  study,  10  of  the  DAS 
were  deployed.  1000  samples  at  1  msec  sample  interval  were  recorded  for  all 
data  sets.  The  DAS  have  a  16  bit  analog  to  digital  converter  with  an  effective 
dynamic  range  of  90  dB.  The  system  has  7  selectable  gains  from  1  to  8192. 
The  geophones  were  Mark  Products,  Inc.  L-410H's.  The  geophones  have  a 
corner  frequency  of  IC*  and  a  damping  of  0.6.(Figure  8).  Sixty  horizontal 
geophones  were  orien  f  '  ansverse  to  the  propagation  path  for  the  SWIG 


le  ai 


mm 

■1 

■ 

■ 

■ 

■ 

1 

1 

1 

■ 

■ 

RESPONSE  CURVE 

L-410  CEOPHOnE 

10.0  HZ  330  OHM  COIL 

■ 

■ 

■ 

1 

1 

1 

II 

1 

mnii 

■I 

■ 

1 

o 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

i 

1 

1 

g 

■ 

■ 

_ 1 _ 

1  1  t 

1 

•  2  • 

t  t 

1  1 

1 

mr 

Z] 

1“ 

u_ 

1 

5r 

1  1 

r' 

1  1 

I _ i _ 

r 

L_ 

“ 

“ 

I  . 

•  (/) 

H 

L_i_ 

1 

L 

■ 

■ 

n 

wu 

SSSssi* 

1 

as 

1 

i 

1! 

\wmmamam 

.^l 

■1 

■ 

■ 

■ 

■ 

1 

m> 

m 

1 

i 

1 

1 

1 

J _ ( 

b 

CL 

.  ►-  . 

i 

■ 

■ 

1 

1 

1 

n 

■■ 

■ 

■ 

3 

O 

1 

■ 

■ 

■n 

1 

1  i 

P 

UHVE  SHUNT  DAMPING 

A  OPEN  0  475 

B  23430  OHMS  0.500 

C  4422  OHMS  0  SCO 

D  2310  OHMS  0.7C0 

1 

1 

1 

I 

1 

m 

1 

1 

FRJ 

=QUE 

^■1 

NCY- 

— HEnTZ 

r 

L_ 

D^«n<0a90  O  OOOOOO  o  cc 

^  CM  n^mtocoo  o  oo 

^  CM  n  ^ 


Figure  8.  Instrument  response  of  L-410  geophones  used  in  this 
study. 
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source  and  sixty  phones  were  oriented  vertically  tor  the  Betsy  source.  Figure  9 
illustrates  the  acquisition  layout.  Both  the  SWIG  (Rgure  6)  and  Betsy  (Figure  7) 
data  sets  were  acquired  with  the  geophones  spaced  at  one  halt  meter  incre¬ 
ments  and  covering  the  range  from  5  to  34.5  meters.  The  sources  were  moved 
by  30  meters  so  that  source-receiver  coverage  was  extended  to  35-64.5  meters. 
This  data  Is  plotted  as  Figure  10  tor  the  SWIG  source  and  Figure  11  tor  the 
Betsy.  The  sources  were  moved  one  more  time  to  extend  the  coverage  trom  65 
to  94.5  meters.  The  tinal  data  sets  are  plotted  in  Figure  12  tor  the  SWIG  source 
and  Figure  13  tor  the  Betsy  source.  The  acquisition  geometry  ot  the  6  data  sets 
are  summarized  in  Table  2.  The  different  arrays  ot  data  acquired  in  this  study 
are  identitied  by  the  spacing  between  the  source  and  the  nearest  receiver. 

The  spectral  data  at  5  meter  increments  tor  the  SWIG  source  are  repro¬ 
duced  in  Figure  14.  These  spectra  with  a  spectral  peak  near  25  Hz  and  data 
above  the  background  noise  to  beyond  100  Hz  demonstrate  the  broad  band  na¬ 
ture  ot  the  source.  The  low  Q  materials  in  the  tormations  cause  the  higher  tre- 
quencies  to  lose  amplitude  tairly  rapidly.  The  observed  spectra  are  not  smooth, 
but  show  notches  which  are  apparently  a  result  ot  propagation  path  effects.  The 
first  hole  in  the  spectra  occurs  at  68  Hz  in  the  5  meter  data  and  then  decreases 
to  57  Hz  at  10.5  meters,  40  Hz  at  15  meters  and  even  lower  with  increasing 
range.  The  decreasing  frequency  ot  this  Interference  hole  with  increasing  range 
may  represent  an  increasing  time  separaJon  of  two  interfering  signals  with 
range. 
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Rgure  9.  Acquisition  geometry  for  both  SWIG  and  Betsy  dats  sets. 
Sources  were  located  at  3  positions:  5, 35  and  65  meters  from  the 
near  geophone.  Note  the  distances  from  the  far  geophones.  The 
far  phone  for  the  5  meter  data  set  is  34.5  meters,  almost  identical 
to  near  phone  distance  for  next  source  at  35.0  meters  although  the 
geology  sampled  by  the  paths  is  quite  different. 
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Figure  12.  Swig  data  with  source  offset  at  65.0  m  from  near  receiver. 
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Table  2 

ACQUISITION  GEOMETRY  OF  6  DATA  SETS 
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Data  set  name 

Source 

Receivers 

5  meter  SWIG 

SWIG  located  5  meters 
from  nearest  receiver 

60  transverse  oriented 
geophones  in  linear 
array  spaced  at  0.5 
meter  increments  over 
interval  5.0  to  34.5 
meters. 

35  meter  SWIG 

SWIG  located  35  meters 
from  nearest  receiver 

60  transverse  oriented 
geophones  in  linear 
array  spaced  at  0.5 
meter  increments  over 
interval  35.0  to  64.5 
meters. 

65  meter  SWIG 

SWIG  located  65  meters 
from  nearest  receiver 

60  transverse  oriented 
geophones  in  linear 
array  spaced  at  0.5 
meter  increments  over 
interval  65.0  to  94.5 
meters. 

5  meter  Betsy 

Betsy  located  5  meters 
from  nearest  receiver 

60  vertical  oriented 
geophones  in  linear 
array  spaced  at  0.5 
meter  increments  over 
interval  5.0  to  34.5 
meters. 

35  meter  Betsy 

Betsy  located  35  meters 
from  nearest  receiver 

60  vertical  oriented 
geophones  in  linear 
array  spaced  at  0.5 
meter  increments  over 
interval  35.0  to  64.5 
meters. 

65  meter  Betsy 

Betsy  located  65  meters 
from  nearest  receiver 

60  vertical  oriented 
geophones  in  linear 
array  spaced  at  0.5 
meter  increments  over 
interval  65.0  to  94.5 
meters. 
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Rgure  14.  Spectral  Response  of  single  receiver  SWIG  data  from  5 
to  30  meters.  Clockwise  from  upper  left,  5  meters,  10  meters,  20 
meters,  30  meters,  25  meters,  15  meters. 
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CHAPTER  3 

DATA  PROCESSING  FOR  SURFACE  WAVE  DISPERSION 


The  data  processing  was  done  on  a  Sun  3/60  workstation  in  the  SMU 
geophysics  lab.  The  processing  primarily  utilized  a  system  of  programs  devel¬ 
oped  as  an  ongoing  project  by  Robert  Herrmann  of  St.  Louis  University.  This 
package  contains  interrelated  programs  which  allowed  the  output  from  one  pro¬ 
gram  to  be  used  as  input  to  another.  The  software  descriptions  throughout  this 
section  will  describe  the  Herrmann  implementations  of  the  various  techniques. 

A.  Pre-processing  of  field  data 

The  field  data  was  copied  from  IBM  PC  compatible  disks  to  the  Sun  work¬ 
station  as  ASCII  files.  The  software  processes  data  starting  from  a  frequency 
domain  representation  of  the  seismic  traces.  The  data  was  converted  to  this 
format  while  simultaneously  detrending  ,  or  removing  the  DC  bias.  The  detrend¬ 
ing  was  accomplished  by  taking  the  average  value  of  the  entire  trace  and  sub¬ 
tracting  that  average  value  from  every  sample  of  the  trace.  The  decision  was 
made  to  process  the  seismic  data  as  six  separate  data  sets  of  sixty  traces  each. 
Each  data  set  contained  a  separate  source  firing.  It  appeared  that  there  were 
some  slight  differences  between  the  data  sets  each  time  the  source  was  moved. 
It  also  seemed  reasonable  that  sixty  traces  were  as  many  traces  as  would  be 
needed  for  any  array  type  processing.  The  previous  section  showed  plots  of  the 
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data  after  detrending  for  the  six  data  sets.  One  bad  trace  exists  on  each  data 
set  due  to  a  connection  problem  during  the  acquisition. 

The  initial  processing  and  previous  results  indicated  that  low  frequency  or 
long  period  data  would  provide  the  best  resolution  of  the  deep  structure.  Since 
the  geophones  used  had  a  corner  frequency  of  10  Hz  an  instrument  response 
correciion  was  applied  to  the  data  in  order  to  expand  the  usable  bandwidth  of 
the  data.  A  comparison  of  the  spectra  at  select  distances  before  and  after  in¬ 
strument  correction  (Figure  15)  show  that  the  only  difference  was  a  boosting  of 
the  amplitude  below  about  10  Hz.  The  instrument  corrected  data  is  reproduced 
in  Figures  16  to  21.  The  seismic  records  after  correction  for  instrument  re¬ 
sponse  were  then  processed  to  determine  dispersion  through  various  tech¬ 
niques.  The  most  fundamental  technique  for  determining  group  velocity  disper¬ 
sion  utilized  in  this  study  was  the  multiple  filter  technique. 


1.  Experimental  Technique 

The  Multiple  RIter  Analysis  (MFA)  technique  used  throughout  this  project 
to  derive  group  velocity  dispersion  was  first  proposed  by  Dziewonski  et.  al. 
(1969).  The  group  velocity  is  determined  by  performing  a  narrow  band  filter  on 
the  data  and  searching  for  the  maximum  amplitude  of  the  envelope  function. 
This  procedure  is  repeated  for  each  frequency  over  the  band  of  interest.  The 
group  velocity  can  be  related  to  the  time  for  the  peak  of  the  envelope  by  equa¬ 
tion  (1.7).  Some  advantages  of  the  MFA  over  the  moving  window  analysis  are 
documented  by  Dziewonski  et  ai.  (1969).  These  advantages  include  better  fre¬ 
quency  resolution,  no  intrinsic  period  shift  and  less  computer  processing  time. 
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Figure  15.  Comparison  of  amplitude  spectrums  of  SWIG  data. 
Data  on  the  left  is  the  original  data.  Data  on  the  right  is  with 
instrument  correction.  From  top  to  bottom,  distance  is  65  meters, 
75  meters  and  85  meters. 
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Figure  18.  SWIG  data  after  Instrument  response  correction  for  source 
offset  of  65  m  from  near  receiver. 
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Figure  20.  Betsy  data  after  instrument  response  correction  for  source 
offset  of  35  m  from  near  receiver. 
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Figure  21 .  Betsy  data  after  instrument  response  correction  for  source 
offset  of  65  m  from  near  receiver. 
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Filtering  is  accomplished  with  a  Gaussian  filter  which  is  set  to  zero  outside 
a  given  bandwidth  and  is  tapered  inside  this  bandwidth  from  the  center  fre¬ 
quency.  The  range  of  center  frequencies  are  selected  by  the  user  with  the 
bandwidth  controlled  by  one  parameter  a(which  will  be  used  interchangeably 
with  AO).  The  bandwidth  is: 

(3.1)  BAND  =  (3i/a)1/2 

The  default  value  for  a  is  50.27  which  yields  a  value  for  BAND  of  25%. 
This  means  that  for  a  center  frequency  of  100  Hz,  the  non-zero  frequency  range 
would  be  75  to  125  Hz.  The  selected  value  for  a  is  a  compromise  between 
resolution  in  the  frequency  domain  and  resolution  in  the  time  domain.  Some 
simple  tests  of  impulse  response  were  made  to  quantify  this  compromise.  An 
impulse  data  set  was  created  which  contained  one  thousand  samples  at  1  msec 
sample  rate.  The  value  of  all  samples  was  set  to  zero  except  for  one  non-zero 
sample  at  0.5  seconds.  The  trace  was  then  filtered  from  1  to  50  Hz  using  a 
Gaussian  bandpass  filter  with  various  values  of  a.  Figure  22.  shows  the  impulse 
response  in  the  time  domain  versus  frequency  for  a  equal  to  5  or  79%  bandwidth 
in  the  frequency  domain.  The  trace  at  approximately  75  Hz  is  the  input  impulse 
trace.  The  corresponding  envelope  appears  at  the  bottom  of  the  page.  Figure 
23  contains  the  same  plots  for  AO  set  to  75  for  comparison. 

Time  domain  sidelobes  from  the  bandwidth  increase  as  frequency  de¬ 
creases  or  the  period  increases.  Below  10  Hz  the  sidelobes  become  so  large  it 
is  difficult  to  discriminate  between  closely  spaced  events  in  time.  The  envelope 
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Figure  22.  Impulse  response  in  time  domain  of  bandpass  filters 
for  different  center  frequencies,  a  or  AO  equals  5  in  this  case  so 
bandwidth  is  79%.  Trace  at  approximately  75  Hz  is  input  impulse 
trace.  Bottom  plot  Is  envelope  function  of  top  plot. 
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Figure  23.  Impulse  response  in  time  domain  of  bandpass  filters 
for  different  center  frequencies,  a  or  AO  equals  75  in  this  case  so 
bandwidth  is  20%.  Trace  at  approximately  75  Hz  is  input  impulse 
trace.  Bottom  plot  is  envelope  function  of  top  plot. 
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function  is  a  Gaussian  shaped  curve  in  the  time  domain,  representing  one 
smooth  output  event  for  one  input  event. 

Figure  24  displays  the  frequency  domain  response  of  the  Gaussian  filter 
for  various  values  of  a  for  a  center  frequency  of  50  Hz.  For  a  value  of  a  equal  to 
75  or  a  bandwidth  of  20%,  it  can  be  seen  that  the  filter  is  very  tight  in  the  fre¬ 
quency  domain  yielding  a  good  frequency  resolution.  Unfortunately,  the  side- 
lobes  are  so  large  in  the  time  domain  that  it  becomes  difficult  to  separate  differ¬ 
ent  velocity  events,  especially  at  low  frequencies. 

Resolution  in  the  time  domain  or  separation  in  time  between  two  events  is 
related  to  the  bandwidth  in  the  frequency  domain  and  the  period  T  (Herrmann, 
1988): 

(3.2)  Time  Resolution  =  T/Band. 

Table  3  summarizes  the  results  of  the  impulse  response  tests,  and  dem¬ 
onstrates  the  tradeoff  between  frequency  resolution  and  time  domain  resolution. 
This  information,  along  with  processing  panels  of  data  with  various  values  of  AO, 
was  used  to  select  AO  for  the  rest  of  the  project.  A  value  of  25  was  used  for  all 
of  tho  data  in  this  project.  In  theory  this  should  only  allow  the  determination  of 
group  velocity  when  interfering  events  are  separated  by  100  msec  or  more. 

The  bandwidth  of  the  acquired  data  used  in  this  study  is  approximately  10 
to  30  Hz.  Some  effort  was  expended  to  recover  quality  data  below  10  Hz,  in 
eluding  correction  for  instrument  response.  It  appears  that  this  was  not  effica¬ 
cious  since  the  major  events  were  separated  by  less  than  200  msec  at  this  dis¬ 
tance.  In  order  for  this  technique  to  work  at  these  frequencies,  it  appears  that 
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Figure  24.  The  impulse  response  in  frequency  domain  at  center 
frequency  of  50  Hz  for  various  values  of  AO.  Clockwise  from  upper 
left,  AO  is  5. 10.  35,  75,  50  and  25. 
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data  would  have  to  be  acquired  at  much  further  source-receiver  offsets  which 
would  provide  larger  time  separation  between  arrivals. 

2.  Experimental  Results 

The  MFT  was  applied  to  all  six  data  sets.  Displays  of  the  MFT  analysis  of 
the  SWIG  generated  Love  wave  data  at  approximately  5  m  intervals  are  pre¬ 
sented  in  Figures  25  to  45.  The  bottom  plot  contains  the  amplitude  contours 
across  the  various  frequencies  on  the  x-axis.  The  y-axis  is  time  scaled  by  group 
velocity  according  to  equation  (1.7).  The  input  trace  is  presented  to  the  far  right 
scaled  in  linear  time  and  also  in  linear  velocity  for  comparison  to  the  contoured 
results.  The  four  largest  peaks  of  the  envelope  function  at  each  frequency  are 
noted  and  marked  on  the  plot  with  the  largest  amplitude  marked  with  a  square, 
second  largest  with  a  circle,  third  with  a  triangle  and  the  fourth  peak  with  a  plus 
sign.  These  amplitudes  are  scaled  relatively  and  presented  in  the  frequency 
domain  in  the  top  plot. 

Some  general  things  can  be  noted  by  studying  figure  25.  The  trace  at  the 
far  right  shows  that  there  is  a  small  period  of  time  in  which  most  of  the  energy 
arrives.  The  second  trace  shows  that  the  peak  amplitude  of  the  trace  arrives 
with  a  velocity  of  approximately  90  m/s.  A  dispersion  curve  could  be  drawn  con¬ 
necting  the  squares  from  3  Hz  to  60  Hz.  The  dispersion  does  not  match  the 
classical  shape  we  would  expect  of  velocity  decreasing  with  frequency.  This  dif¬ 
ference  is  due  in  part  to  the  fact  that  the  body  waves  and  surface  waves  have 
not  had  enough  time  to  separate  at  this  range.  Another  reason  could  be  that  the 
Love  waves  have  not  had  sufficient  time  to  stabilize.  One  other  feature  that  can 
be  seen  in  the  spectral  plot  at  the  top  is  that  the  amplitude  of  the  main  event  is 


51 


over  an  order  of  magnitude  larger  than  the  next  largest  event  over  a  broad  band. 
This  implies  that  the  data  set  contains  information  from  8-60  Hz. 

Figure  26  shows  the  bandwidth  dropping  to  8-55  Hz  at  a  source-receiver 
spacing  of  10.5  m.  There  is  no  indication  of  dispersion  with  the  main  body  of 
energy  arriving  at  100  m/s.  At  20  ms  (Figure  28)  there  is  some  indication  of  dis¬ 
persion.  Group  velocity  for  the  event  of  interest  is  about  150  m/s  at  13  Hz  and 
drops  to  about  100  m/s  at  45  Hz.  The  spectral  plot  shows  two  strong  events  that 
cross  at  60  Hz.  This  information  can  be  used  to  help  determine  continuity  of  an 
event  on  the  velocity  plot.  In  this  case,  the  main  event  would  be  a  square.  The 
dispersion  curve  would  be  drawn  connecting  squares  from  8  Hz.  to  about  60  Hz. 
The  dispersion  curve  would  be  drawn  connecting  circles  from  that  point  forward. 
This  technique  would  be  most  useful  if  the  velocities  of  two  events  converge  for 
a  small  interval  and  then  diverge. 

At  30  m  (Figure  30)  the  trace  can  be  seen  to  clearly  contain  two  events. 
The  first  event  of  velocities  180  to  300  m/s  is  clearly  the  direct  and  refracted 
shear  body  waves.  The  largest  amplitude  event  on  the  trace  arrives  directly  af¬ 
terwards,  corresponding  to  a  velocities  between  100  and  180  m/s.  The  event 
has  several  cycles  with  the  longest  period  portion  arriving  earliest  and  is  inter¬ 
preted  as  a  Love  wave  with  normal  dispersion.  As  the  source-receiver  spacing 
increases,  the  Love  wave  portion  generally  becomes  more  widely  separated  in 
time  from  the  body  waves  or  direct  arrivals.  The  Love  wave  portion  of  the  seis¬ 
mic  trace  also  shows  increased  dispersion  with  offset. 

The  group  velocity  determined  from  the  seismic  trace  is  an  average  over 
the  spatial  interval  between  the  source  and  receiver.  This  implies  that  for  any 
given  receiver,  say  at  80  m,  the  next  receiver  should  have  very  similar  results.  It 
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can  be  seen  from  analyzing  Figure  37  and  38  that  the  results  at  80  and  85  m  are 
very  similar,  in  spite  of  a  being  relatively  noisy.  If  the  site  were  horizontally 
stratified  and  homogeneous,  two  receiver  locations  0.5  m  apart  should  be  ex¬ 
tremely  similar. 

There  exist  two  portions  of  the  data  set  that  are  at  approximately  the 
same  source  receiver  spacing,  but  cover  a  different  interval  of  ground.  When 
the  data  was  acquired,  the  source  was  moved  back  30  m  at  a  time  with  the  re¬ 
ceivers  left  in  the  same  60  locations.  Figure  9  demonstrates  that  the  source  at 
35  m  from  the  near  receiver  and  the  source  at  34.5  m  from  the  far  receiver  cover 
a  similar  distance,  although  the  travel  path  for  these  two  stations  overlaps  by 
only  5  m.  Figure  31  and  32  are  analyzed  at  34.5  and  35  m,  and  exhibit  different 
Love  wave  dispersion  as  well  as  different  amplitude  levels  for  the  entire  seismic 
trace.  Differences  between  these  two  traces  is  most  pronounced  at  the  higher 
frequencies.  The  two  data  sets  at  64.5  and  65  m  (Figures  38  and  39) 
demonstrate  even  larger  differences.  The  data  set  at  64.5  m  has  almost  4  full 
cycles  apparent  on  the  raw  trace  of  the  Love  wave.  The  Love  wave  in  the  trace 
at  65  m  consists  of  only  about  two  cycles  .  This  alone  is  enough  to  indicate  that 
the  two  intervals  imaged  by  the  Love  wave  are  geologically  different.  The  dis¬ 
persion  results  are  very  similar  between  10  and  30  Hz.  The  high  frequency 
components  show  the  most  variability  over  the  entire  data  set. 

Figure  46  is  a  composite  diagram  of  the  group  velocity  dispersion  for  the 
Love  wave  on  the  5  meter  data  set  at  5  m  increments.  The  velocity  increases  for 
each  dispersion  curve  as  the  range  increases  from  5  to  34.5  m.  This  result  is 
consistent  with  Mari's  (1984)  statement  that  it  would  take  an  offset  distance  of  4 
to  6  times  the  deepest  layer  imaged  for  the  Love  wave  to  stabilize.  The  velocity 


Rgure  26.  MFT  analysis  at  10.5  m  for  SWIG  data 


Figure  31 .  MFT  analysis  at  34.5  m  tor  SWIG  data 
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Figure  34.  MFT  analysis  at  45.0  m  for  SWIG  data 
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Figure  35.  MFT  analysis  at  50.0  m  for  SWIG  data 
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Figure  36.  MFT  analysis  at  55.0  m  for  SWIG  data 
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Figure  42.  MFT  analysis  at  80.0  m  for  SWIG  data 
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Figure  46.  Interpreted  dispersion  curves  for  Love  wave  group 
velocity  from  multiple  filter  analysis  at  range  of  5  to  34.5  m.  Source 
at  5  m  from  near  receiver. 
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Figure  48.  Interpreted  dispersion  curves  for  Love  wave  group 
velocity  from  multiple  filter  analysis  at  range  of  65  to  94.5  m. 
Source  at  65  m  from  near  receiver. 
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still  shows  some  Increase  over  the  Interval  Irom  35  to  64.5  m  (Figure  47),  but 
substantially  less  change  than  observed  In  the  5  to  34.5  m  data  set.  This  result 
Implies  that  the  Love  wave  Is  not  fully  stabilized  by  35  m,  but  Is  nearly  so  at  65 
m.  The  final  composite  plot,  Figure  48,  shows  the  MFA  Interpreted  dispersion 
curves  from  65  to  94.5  m.  The  values  are  almost  coincident  between  10  and  25 
Hz.  The  Love  wave  appears  stabilized  and  usable  for  making  geological  inter¬ 
pretations  at  these  offsets. 

The  Betsy  data  sets  covers  the  same  intervals,  but  are  slightly  more 
complicated  since  they  contains  both  a  fundamental  mode  as  well  as  the  first 
higher  mode  Rayleigh  wave.  In  addition,  the  temporal  separation  of  the  events 
is  smaller,  making  analysis  more  difficult.  The  multiple  filter  processing  of  the 
Betsy  data  set  is  contained  in  the  Appendix. 

C.  Array  Processing  of  Group  Velocity 
1.  Experimental  Technique 

The  Herrmann  (1987)  software  package  used  in  this  study  implements  a 
technique  from  Barker  (1988)  which  makes  use  of  an  array  of  seismograms  to 
determine  group  velocity.  This  technique  uses  the  multiple  filter  algorithms  pre¬ 
viously  documented.  Each  trace  is  bandpassed  with  a  series  of  Gaussian  filters. 
The  envelope  of  each  resulting  bandpassed  trace  is  then  computed.  The  addi¬ 
tional  step  involved  in  array  processing  is  to  stack  the  normalized  envelopes 
along  lines  of  constant  velocities.  As  Herrmann  notes,  this  step  is  not  a  slant 
stack,  since  only  a  single  point  results  from  stacking  along  each  velocity.  Up  to 
four  peaks  of  this  stacked  envelope  function  are  selected  to  yield  group  veloci¬ 
ties  of  up  to  four  modes  at  a  particular  frequency.  This  processing  is  then  fol- 


lowed  by  a  more  detailed  analysis,  in  which  each  original  envelope  trace  is 
searched  around  each  of  the  four  average  velocities  to  obtain  the  individual 
group  velocities,  which  generally  differ  from  the  initial  estimate.  The  individual 
velocities  are  used  to  compute  a  mean  group  velocity  and  standard  error.  These 
velocities  represent  an  average  response  of  the  material  being  imaged  by  the 
different  source  receiver  pairs.  As  with  the  multiple  filter  analysis,  the  different 
modes  must  be  separated  in  time  by: 

(3.3)  T{A0/3t)1/2, 

where  T  is  the  period  of  the  current  bandpass  filter.  If  several  modes  occur  too 
near  one  another,  the  results  of  this  analysis  would  indicate  only  one  mode  pre¬ 
sent.  The  group  velocities  would  thus  be  in  error  due  to  interference  phenom¬ 
ena. 

One  additional  feature  of  this  array  analysis  is  that  it  can  be  used  to  calcu¬ 
late  the  anelastic  attenuation  coefficient,  y.  The  peak  amplitudes  of  the  unnor¬ 
malized  envelopes  are  input  into  a  regression  analysis  of  logarithmic  amplitude 
versus  offset  which  assumes  a  model  of  the  form: 

(3.4)  In  A(r,a))  ss  In  Aq  -  7R(to)r 

The  option  exists  for  making  a  geometrical  spreading  correction  before  calculat¬ 
ing  the  attenuation  coefficients.  The  anelastic  coefficients  determined  here  can 
be  used  to  invert  for  Q.  This  procedure  will  be  used  later  in  the  study  to  limit  the 
uncertainty  in  the  inversion  of  the  velocity  dispersion. 
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2.  Experimental  Results 

The  Barker  technique  tor  group  velocity  determination  was  applied  to  the 
SWIG  data  in  3  groups.  The  value  for  AO  used  for  this  work  was  25,  correspond¬ 
ing  to  a  bandwidth  of  35%,  the  same  value  that  was  used  for  the  multiple  filter 
analysis.  The  eleventh  trace  on  all  three  sections  of  the  SWIG  data  set  was 
noisy  and  was  not  included  in  the  processing.  One  of  the  useful  features  of  this 
type  of  processing  is  that  the  stacking  is  across  equal  velocities.  Uneven  spatial 
sampling  can  thus  be  accommodated. 

The  outputs  of  this  processing  on  ti.e  three  sections  of  the  SWIG  data 
(Rgures  49  through  51)  contain  contours  of  the  amplitude  corresponding  to  99%, 
90%,  80%,  70%,  60%  and  50%  of  the  maximum  amplitude  of  the  envelopes.  As 
with  the  MFT,  the  largest  amplitude  peak  is  marked  by  squares,  second  largest 
by  circles,  third  largest  by  triangles  and  the  fourth  largest  by  plus  signs.  The 
largest  event  over  the  interval  10  to  30  Hz  corresponds  to  the  fundamental  mode 
Love  wave  with  its  frequency  content  decreasing  with  distance.  The  near  data 
set  appears  coherent  from  6  Hz  to  80  Hz  with  velocities  ranging  from  about  95  to 
150  m/s.  The  group  velocity  shows  more  dispersion  and  higher  velocities  with 
distance.  This  effect  appears  to  be  the  result  of  some  near  field  effects  com¬ 
bined  with  the  fact  that  it  takes  some  time  (or  distance)  for  the  surface  waves  to 
be  set  up,  especially  at  the  longer  periods.  The  velocity  dispersion  of  the  3  data 
sets  (Figure  52)  show  no  overlap  even  at  the  higher  frequencies.  This  difference 
certainly  implies  that  an  offset  greater  than  65  m  would  be  necessary  for  consis¬ 
tent  results.  Figures  53  to  59  are  overlays  of  the  single  station  MFT  results  with 
the  array  results  for  the  data  set  from  65  to  94.5  m.  The  results  are 
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Figure  49.  Barker  technique  airay  processing  of  the  SWIG  data 
set  for  group  velocity  over  the  interval  from  5  to  34.5  m 


0.00 


13.33 


66.67 


26.67 


40.00 


53.33 


FREQUENCY  (HZ) 


Rgure  51 .  Barker  technique  array  processing  of  the  SWIG  data 
set  for  group  velocity  over  the  interval  from  65  to  94.5  m 
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Rgure  52.  Overlay  of  the  Barker  technique  group  velocity 
dispersion  from  5, 35  and  65  meter  SWIG  data  sets 
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Figure  53.  Single  station  processing  at  65.0  m  (triangles)  versus 
array  processing  using  Barker  technique(circles)  of  65  meter  SWIG 
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Figure  54.  Single  station  processing  at  70.5  m  (triangles)  versus 
array  processing  using  Barker  technique(circles)  of  65  meter  SWIG 
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very  consistent  up  to  29  Hz.  The  results  above  29  Hz  show  an  increased 
amount  of  scatter. 

Examination  of  the  endpoints  of  the  data  where  the  sources  were  moved 
illustrate  some  significant  aspects  of  the  structure.  Figure  60  shows  the  single 
station  analysis  of  group  velocity  at  34.5  and  35.0  m.  The  single  station  results 
are  also  contrasted  to  the  corresponding  array  velocities  at  5  to  34.5  m  and  35  to 
64.5  m.  Although  there  is  only  0.5  meter  difference  in  range  between  the  two 
single  station  measurements,  there  is  a  substantial  difference  in  the  group  veloc¬ 
ity  curves.  The  two  velocities  are  almost  the  same  from  6  to  12  Hz  and  differ  by 
about  10  percent  above  20  Hz.  The  acquisition  involved  leaving  the  geophones 
stationary  and  moving  the  sources  to  increase  the  range.  The  group  velocity  is 
the  result  of  the  average  velocity  structure  between  the  source  and  receiver.  In 
this  case,  although  there  is  only  one  half  meter  change  In  range,  there  Is  only  5.0 
m  overlap  in  the  interval  covered.  This  leads  to  the  conclusion  that  the  differ¬ 
ence  between  the  two  group  velocities  is  due  to  lateral  changes  in  the  velocity 
field.  The  fact  that  these  velocities  are  different  at  higher  frequencies  and  the 
same  at  low  frequencies  indicates  that  the  velocity  structure  is  different  for  the 
shallower  layers,  but  the  same  for  the  deeper  structure. 

These  conclusions  can  be  checked  by  comparing  the  results  at  64.5  m 
and  65.0  m  in  Figure  61 .  Once  again,  the  two  ranges  are  0.5  m  apart,  but  the 
covered  interval  differs  by  30  m.  The  two  dispersion  results  are  similar  at  low 
frequencies,  but  differ  by  25  percent  at  frequencies  from  20  to  40  Hz.  The  other 
interesting  result  is  that  both  of  these  single  station  estimates  are  comparable  to 
the  corresponding  array  results.  The  single  station  dispersion  results  at  64.5  m 
is  only  slightly  higher  than  the  array  or  average  results  from  35  to  64.5  m.  The 
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single  station  group  velocity  at  65  m  matches  almost  perfectly  the  array  average 
over  the  interval  65  to  94.5  m. 

As  noted  earlier,  the  array  processing  provides  the  opportunity  to  estimate 
attenuation  coefficients.  These  estimates  Increase  with  frequency  (Figure  62)  al¬ 
though  there  is  some  variability  in  the  curves.  The  three  data  sets  show  similar 
attenuation  at  the  lower  frequencies,  between  8  and  20  Hz.  As  with  all  the  data, 
there  is  greater  scatter  at  frequencies  above  30  Hz. 

The  variability  of  the  velocity  and  attenuation  estimates  at  high  frequen¬ 
cies  is  a  process  that  needs  to  be  understood  If  one  is  to  fully  utilize  the  results. 
One  probable  explanation  for  this  variation  is  that  the  assumption  of  horizontally 
stratified  layers  is  being  violated  in  the  shallowest  part  of  the  subsurface.  At  a 
frequency  of  30  Hz,  with  the  shallowest  layer  shear  velocity  of  150  m/s,  one 
wavelength  is  about  5  m.  If  the  shallowest  layer  was  varying  smoothly  from  1 .5 
to  4  m  in  thickness,  the  dispersion  results  for  the  higher  frequencies  could  show 
substantial  variance.  Bogaards  (1989)  indicated  that  the  high  frequency  vari¬ 
ations  could  be  the  result  of  scattering  off  of  caliche  lenses. 

The  high  frequency  variability  indicates  more  geologic  complexity  in  the 
shallowest  environments  than  was  expected  at  the  start  of  this  study.  The  con¬ 
sistencies  at  lower  frequencies,  or  longer  wavelengths,  indicates  that  the  deeper 
horizons  are  more  geologically  uniform  than  the  shallowest  environment. 

The  body  wave  portion  of  the  wavefield  in  the  array  processing  shows 
some  variance  in  velocity  versus  frequency,  but  is  centered  around  200  m/s  for 
the  nearest  data  set.  The  body  wave  velocity  increases  with  offset.  This  result 
is  consistent  with  a  normal  refraction  survey  with  an  increase  in  velocity  with 
depth  for  various  layers. 
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Figure  60.  Comparison  of  Multiple  filter  analysis  and  array  techniques  at  34.5  m 
and  35.0  m.  34.5  m  single  station  processing  represented  by  X;  35.0  m  single 
station  by  circle,  5  to  34.5  m  array  by  +,  and  35  to  64.5  m  array  analysis  by 
circle; 


GAMMA  (1/M) 


95 


0.240 


0.200 


0. 160 


0.120 


0.080 


0.040 


0.000 


SWIG  data  from  5,35  &  65  meters 


FREQUENCY  (HZ) 


Figure  62.  Attenuation  coefficients  from  array  processing(Barker 
Technique)  of  all  three  SWIG  source  positions. 
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The  Betsy  source  data  set  is  more  complicated,  but  the  basic  results  are 
similar.  The  Rayleigh  data  includes  a  fundamental  mode  and  one  higher  mode. 
Both  modes  show  more  clearly  defined  dispersion  and  higher  group  velocities  as 
the  source  is  moved  further  away.  This  relationship  is  consistent  with  the  earlier 
observation  that  some  time  or  distance  had  to  elapse  before  the  surface  waves 
are  fully  developed.  The  65  meter  data  set  has  a  further  complication  in  that  the 
air  blast  at  300  m/s  couples  to  the  higher  mode  data  at  high  frequencies. 

Figures  63  to  65  show  the  Barker  technique  results  for  5,  35  and  65  meter 
Betsy  data  sets.  The  energy  for  all  of  these  data  sets  fall  into  two  distinct  bands. 
The  faster,  higher  frequency  band  is  the  first  higher  mode  Rayleigh  wave.  There 
also  tends  to  be  a  less  dispersive  portion  at  300  m/s.  This  energy  is  generated 
by  the  air  blast  from  the  source.  The  interpreted  results  of  the  three  fundamental 
mode  data  sets  are  shown  in  Figure  66.  The  35  meter  data  set  breaks  up  at 
about  25  Hz  and  shows  some  sort  of  Airy  phase  at  about  18  Hz.  This  analysis  is 
consistent  with  the  results  from  Simila  (1982).  Figure  67  contains  the  first  higher 
mode  dispersion  curves  for  all  3  data  sets.  The  65  meter  data  set  appears  to  be 
strongly  biased  by  the  air  blast  above  45  Hz,  and  possibly  may  at  lower  fre¬ 
quencies  as  well.  The  attenuation  coefficients  for  the  fundamental  modes 
(Rgure  68)  are  fairly  consistent  at  all  distances.  The  higher  mode  attenuation 
coefficients  (figure  69)  show  negative  attenuation  or  an  increase  in  amplitude 
with  distance  for  the  65  meter  data  set.  This  physically  impossible  result  ap¬ 
pears  to  be  due  to  the  coupling  of  the  air  blast  energy  back  into  the  wavefield. 
The  5  and  35  meter  array  data  sets  have  more  normal  appearing  attenuation. 
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Figure  65.  Array  processing  for  group  velocity  of  Betsy  data  set  at 
65  m.  Highest  energy  corresponds  to  squares,  second  highest 
circles,  third  highest  triangles,  and  fourth  highest  plus  signs. 
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Figure  66.  Overlay  of  the  interpreted  group  velocity  of  Rayleigh 
wave  fundamental  mode  from  5,35  and  65  meter  ^tsy  data  sets 
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Figure  67.  Overlay  of  the  interpreted  group  velocity  of  Rayleigh 
wave  first  higher  mode  from  5,35  and  65  meter  Betsy  data  sets 
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Figure  68.  Overlay  of  the  interpreted  attenuation  coefficients  of 
Rayleigh  wave  fundamental  mode  from  5,35  and  65  meter  Betsy 
data  sets 
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Rgure  69.  Overlay  of  the  interpreted  attenuation  coefficients  of 
Rayleigh  wave  first  higher  mode  from  5,35  and  65  meter  Betsy 
data  sets 


D.  Phase  Matched  Filter 
1.  Experimental  Technique 
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A  complementary  technique  for  calculating  velocity  dispersion  is  phase- 
matched  filtering  or  PMF  as  documented  by  Herrin  and  Goforth  (1977)  and  Go¬ 
forth  and  Herrin  (1979).  The  PMF  procedure  starts  with  an  initial  estimate  of 
dispersion,  generally  from  the  multiple  filter  technique.  An  iterative  technique  Is 
used  to  find  and  apply  a  filter  that  is  phase-matched  for  extracting  a  particular 
mode  or  arrival  of  a  surface  wave.  This  phase-matched  filter  is  then  used  to  cal¬ 
culate  the  final  dispersion  curve. 

Cross-correlation  of  a  signal,  s(t),  with  a  filter,  f(t),  can  be  represented  in 
the  frequency  domain  as: 

(3.5)  s(t)  ®  f(t)  =  |S(o>)||F(co)|expi[a(a))  -  (t>(to)], 

where  5(q})  and  <)>(o))  are  the  phase  spectrums  of  s(t)  and  f(t),  respectively.  If  f(t) 
is  chosen  such  that  the  Fourier  phase  is  the  same  as  that  of  s(t),  then  f(t)  is  a 
phase-matched  filter  with  respect  to  s(t).  The  output  of  the  cross-correlation  is 
then  (S(a))|F(a))|  and  will  be  an  even  function  in  the  time  domain  as  is  the  auto¬ 
correlation  function.  This  output  is  called  the  pseudo-autocorrelation  function  or 
PAF.  The  PAF  depends  only  upon  the  selection  of  the  amplitude  spectrum  of 
the  phase  matched  filter  if  the  phase  of  the  filter  and  signal  match.  |F(u))|  is  se¬ 
lected  to  be  equal  to  1 ,  which  means  that  the  PAF  will  tend  to  have  a  minimum 
width  spike  centered  around  the  time  of  the  event.  This  event  can  be  windowed 
in  the  lag  domain  or  correlation  domain  and  an  FFT  applied.  This  process  sepa¬ 
rates  out  the  portion  of  the  data  that  is  of  interest  and  a  more  accurate  phase 
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matched  filter  can  be  calculated.  The  process  is  repeated  until  the  solution  con¬ 
verges  and  the  phase  of  the  required  mode  is  determined.  The  final  output  is 

the  group  velocity  as  a  function  of  period  computed  from  the  phase.  The  group 
velocity  is  calculated  from  the  group  delay  which  is  related  to  the  Fourier 

phase  by: 

(3.6)  T,(“)  =  ^^. 

^  do) 

where  (|»(c))  is  the  matched  phase  output  from  the  PMF  program.  The  advantage 
of  using  the  PMF  is  that  the  final  results  are  less  contaminated  by  higher  modes 
or  multi-pathed  arrivals. 


2.  Experimental  Results 

The  phase  matched  fitter  technique  is  applied  to  a  single  source  receiver 
pair  at  a  time.  Like  any  model  based  technique,  phase  matched  filtering  results 
are  dependent  on  the  quality  of  the  input,  which  is  a  group  velocity  dispersion 
curve.  The  most  consistent  dispersion  curve  inputs  for  the  SWIG  data  set  are 
those  developed  using  the  Barker  technique  on  the  65  meter  data  set.  The 
computed  group  and  phase  velocities  from  the  phase  matched  filter  contrasted 
to  the  array  technique  results  are  shown  in  Figure  70  to  76  at  approximately  5 
meter  increments. 

The  PMF  results  are  consistent  with  the  multiple  filter  results.  The  single 
phone  group  velocity  deviates  from  the  array  averages  above  25  Hz.  Belo.v  25 
Hz  the  results  of  all  three  techniques  are  comparable.  The  PMF  technique 
tended  to  be  unstable  at  the  higher  frequencies.  This  instability  may  be  a  result 
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65m  SWIG  Barker  vs  65.0m  PMF 
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Figure  70.  Comparison  of  array  derived  group  velocity  with  PMF 
results  at  65  m.  PMF  results  are  indicated  by  circles.  The  top 
curve  is  phase  velocity.  Barker  technique  results  are  represented 
by  triangles. 
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Figure  71 .  Comparison  of  array  derived  group  velocity  with  PMF 
results  at  70.5  m.  PMF  results  are  indicated  by  circles.  The  top 
curve  is  phase  velocity.  Barker  technique  results  are  represented 
by  triangles. 
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Rgure  72.  Comparison  of  array  derived  group  velocity  with  PMF 
results  at  75  m.  PMF  results  are  indicated  by  circles.  The  top 
curve  is  phase  velocity.  Barker  technique  results  are  represented 
by  triangles. 
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Rgure  73.  Comparison  of  array  derived  group  velocity  with  PMF 
results  at  80  m.  PMF  results  are  indicated  by  circles.  The  top 
curve  is  phase  velocity.  Barker  technique  results  are  represented 
by  triangles. 
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Figure  74.  Comparison  of  array  derived  group  velocity  with  PMF 
results  at  85  m.  PMF  results  are  indicated  by  circles.  The  top 
curve  is  phase  velocity.  Barker  technique  results  are  represented 
by  triangles. 
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65m  SWIG  Barker  vs  90.0m  PMF 
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Figure  75.  Comparison  of  array  derived  group  velocity  with  PMF 
results  at  90  m.  PMF  results  are  Indicated  by  circles.  The  top 
curve  is  phase  velocity.  Barker  technique  results  are  represented 
by  triangles. 
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65m  SWIG  Barker  vs  94.5m  PMF 
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Figure  76,  Comparison  of  array  derived  group  velocity  with  PMF 
results  at  94.5  m.  PMF  results  are  indicated  by  circles.  The  top 
curve  is  phase  velocity.  Barker  technique  results  are  represented 
by  triangles. 
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of  other,  interfering  events.  The  MFA  plots  (Figures  34  to  40)  also  indicate  that 
there  is  a  large  amplitude  interference  in  the  frequency  Jomain  above  25  Hz  with 
the  main  body  wave  or  direct  and  refracted  shear  arrivals.  The  velocities  of  the 
body  waves  are  substantially  faster  than  the  velocities  of  the  Love  waves,  180 
versus  320  m/s,  but  the  time  separation  of  these  two  events  is  only  about  1 00 
msec.  The  general  trend  of  the  PMF  results  was  to  give  slightly  lower  velocity 
estimates  than  the  array  technique.  This  result  is  probably  due  to  the  fact  that 
the  PMF  technique  is  less  biased  by  the  faster  body  waves  than  the  other  tech¬ 
niques  used  in  this  thesis. 

The  PMF  technique  generates  a  phase  matched  spectrum  from  which  it  is 
possible  to  create  a  waveform  containing  only  the  surface  wave  being  analyzed. 
Figure  77  contains  a  plot  of  the  extracted  Love  waves  using  this  technique.  The 
traces  all  have  the  classic  dispersed  appearance,  with  a  low  frequency  at  the 
smallest  times  and  high  frequencies  arriving  later.  This  technique  was  used  to 
effectively  subtract  the  body  waves.  These  traces  could  easily  be  subtracted 
from  the  original  traces  to  generate  a  body  wave  trace  without  interfering  surface 
waves. 

The  traces  after  PMF  with  only  the  Love  waves  remaining  were  processed 
through  the  multiple  filter  technique.  The  results  (Figures  78  to  84)  show  that  the 
amplitudes  of  all  other  events  have  been  dropped  by  an  order  of  magnitude. 
The  maximum  frequency  for  the  input  dispersion  curve  to  PMF  was  limited  to  55 
Hz.  No  interpretations  or  inferences  should  be  made  for  the  resulting  data 
above  this  frequency.  The  results  look  quite  coherent  above  this  frequency,  but 
the  low  amplitudes  preclude  any  interpretation. 
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Some  attempts  were  made  at  analyzing  the  65  meter  Rayleigh  wave  data 
set  using  the  PMF  technique.  These  attempts  yielded  inconsistent  results.  One 
problem  with  this  analysis  is  that  the  data  set  tended  to  be  more  complex,  in  part 
because  of  the  existence  of  a  fundamental  and  one  higher  mode.  The  body 
wave  portion  of  the  data  set  appeared  to  contain  more  relative  energy  and  arri¬ 
vals  than  the  same  portion  of  the  SWIG  data  set.  The  higher  velocities  mean 
that  the  events  being  separated  were  closer  together  in  time,  thus  potentially 
being  past  the  resolving  power  of  the  PMF  technique.  The  other  problem  is  that 
the  non-dispersive  air  blast  coupled  back  into  the  Rayleigh  wave,  and  thereby 
distorted  the  dispersion  of  surface  wave  energy. 

E.  Array  Processing  For  Phase  Velocity 
1.  Experimental  Technique 

The  phase  velocity  determinations  using  the  array  data  are  based  on  an 
algorithm  suggested  by  McMechan  and  Yedlin  (1981)  and  implemented  by  Herr¬ 
mann  (1987).  This  technique  relies  on  two  steps  starting  with  a  slant  stack  or 
transform  into  the  t-p  domain.  This  step  is  followed  by  a  Id  transform  or  FFT 
along  lines  of  constant  slowness,  converting  the  inteicept  dimension  to  fre¬ 
quency.  Since  slowness(p)  is  the  inverse  of  velocity,  this  effectively  results  in  a 
direct  representation  of  amplitude  on  a  phase  velocity  versus  frequency(o))  rep¬ 
resentation,  and  is  therefore  called  the  p-omega  technique.  The  maximum  am¬ 
plitudes  correspond  to  dispersion  of  various  modes. 

Mokhtar  et  al  (1988)  describe  an  alternate  method  of  mapping  into  the  p- 
0)  domain  without  doing  a  slant  stack.  The  procedure  is  described  in  detail  in 
Herrmann  (1988)  program  documentation  for  the  p  omega  stack  procedure.  The 
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Figure  77.  PMF  extracted  Love  waves  from  65m  SWIG  data 
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Figure  78.  MFT  processing  of  PMF  stripped  Love  waves  at  65  m. 
Uppe-  plot  is  amplitude  spectrums  corresponding  to  lower  plo* 
group  velocity  dispersion.  Largest  amplitudes  are  squares,  second 
are  circles,  third  are  triangles  and  fourth  amplitude  at  a  frequency 
is  indicated  by  a  plus  sign. 
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Figure  79.  MFT  processing  of  PMF  stripped  Love  waves  at  70.5 
m.  Upper  plot  is  amplitude  spectrums  corresponding  to  lower  plot 
group  velocity  dispersion.  Largest  amplitudes  are  squares,  second 
are  circles,  third  are  triangles  and  fourth  amplitude  at  a  frequency 
is  indicated  by  a  plus  sign. 
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Rgure  80.  MFT  processing  of  PMF  stripped  Love  waves  at  75  m. 
I'poer  plot  is  amplitude  spectrums  corresponding  to  lower  plot 
c .  jijp  velocity  dispersion.  Largest  amplitudes  are  squares,  second 
^  circles,  third  are  triangles  and  fourth  amplitude  at  a  frequency 
IS  indicated  by  a  plus  sign. 
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Figure  80.  MFT  processing  of  PMF  stripped  Love  waves  at  75  m. 
Upper  plot  is  amplitude  spectrums  corresponding  to  lower  plot 
group  velocity  dispersion.  Largest  amplitudes  are  squares,  second 
are  circles,  third  are  triangles  and  fourth  amplitude  at  a  frequency 
is  indicated  by  a  plus  sign. 
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Figure  82.  MFT  processing  ot  PMF  stripped  Love  waves  at  85  m. 
Upper  plot  is  amplitude  spectrums  corresponding  to  lower  plot 
group  velocity  dispersion.  Largest  amplitudes  are  squares,  second 
are  circles,  third  are  triangles  and  fourth  amplitude  at  a  frequency 
is  indicated  by  a  plus  sign. 
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Figure  83.  MFT  prcxiessing  of  PMF  stripped  Love  waves  at  90  m. 
Upper  plot  is  amplitude  spectrums  corresponding  to  lower  plot 
group  velocity  dispersion.  Largest  amplitudes  are  squares,  second 
are  circles,  third  are  triangles  and  fourth  amplitude  at  a  frequency 
is  indicated  by  a  plus  sign. 
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Rgure  84.  MFT  processing  of  PMF  stripped  Love  waves  at  94.5 
m.  Upper  plot  is  amplitude  spectrums  corresponding  to  lower  plot 
group  velocity  dispersion.  Largest  amplitudes  are  squares,  second 
are  circles,  third  are  triangles  and  fourth  amplitude  at  a  frequency 
is  indicated  by  a  plus  sign. 
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seismic  trace  from  each  offset,  rj,  is  represented  in  the  Fourier  domain  by  A(o),rj) 
exp  j<j)j(o)).  The  p-w  stack  is  defined  by  the  relation: 


N 

(3.7)  U(p,oj)  =  ]^C(aj)"’A{£o,r,)exp(jq),)exp(jcopr,)^ 

where  C((o)  =  A(a),r‘|)exp(j(j)-j).  Division  by  C(o))  normalizes  the  spectrum  and 
removes  the  source  phase  if  the  input  trace  is  a  single  mode.  The  actual  spec¬ 
trum  is  assumed  to  be  a  sum  of  multiple(M)  modes  of  surface  waves: 


M 

(3.8)  A(w,r|)exp(j(pj)  =  5^S^(o),rj)expj[(poK(o)) -copok(o))rJ, 

k>1 

where  S|<((o,rj)  is  the  amplitude  of  the  kth  mode  and  the  phase  is  separated  into 
distance-independent  (pok(«))>  and  distance  dependent  (Dpoi^(co)rj  components. 
POk  is  the  inverse  of  the  phase  velocity  of  the  kth  mode  for  frequency  co.  If  the 
input  consists  of  a  single  noise-free  surface  wave  mode,  then  the  quantity  LKp.w) 
will  be  a  maximum  at  p  s  pgk. 

2.  Experimental  Results 

The  software  implementation  of  this  technique  determines  a  maxima  for 
U(p,o)).  A  stack  is  created  over  a  user  specified  number  of  phase  velocities.  Up 
to  four  frequencies  for  which  the  amplitude  is  maximum  are  then  computed  for 
each  velocity.  One  difficulty  that  occurred  was  that  If  the  number  of  stacking  ve¬ 
locities  was  too  large,  there  would  be  too  many  phase  velocities  for  each  fre¬ 
quency.  Figure  85  is  the  result  of  p  omega  processing  of  the  SWIG  data  at  65  m 
stacking  over  100  velocities.  Figure  86  is  the  same  data  set  using  only  40 
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Figure  86.  P-omega  stacking  of  65  meter  SWIG  data.  Nray  or 
number  of  velocities  computed  is  set  to  40. 
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velocities.  Note  that  the  energy  appears  to  die  out  or  become  incoherent  be¬ 
yond  25  Hz.  Both  velocity  stacks  resulted  in  similar  amplitude  contours  in  the  5 
to  25  Hz  range.  The  40  velocity  results  were  easier  to  interpret.  Trace  number 
twelve  was  bad  in  each  data  set  and  so  the  first  twelve  traces  were  excluded 
from  the  p-omega  processing.  The  remaining  48  traces  appeared  to  yield  stable 
results. 

Figures  87  and  88  illustrate  the  40  velocity  processing  of  the  SWIG  data 
at  5  and  35  m.  Both  of  these  data  sets  contain  information  at  higher  frequencies 
than  the  65  meter  data  set.  The  interpreted  phase  velocity  dispersion  curves  are 
overlain  for  comparison  in  Figure  89.  The  most  notable  result  is  that  the  phase 
velocities  are  consistent  over  the  entire  range  of  the  experiment. 

The  p-omega  processing  was  also  applied  to  the  Betsy  data.  The  last 
eighteen  traces  were  not  used  because  of  the  bad  trace.  The  results  of  process¬ 
ing  the  5  meter  data  set  is  presented  in  Figure  90.  As  with  the  array  processing 
for  group  velocity,  the  contours  represent  99  percent,  90  percent,  80  percent,  70 
percent,  60  percent  and  50  percent  of  the  maximum  amplitude.  Two  bands  of 
data  appear  on  the  processing.  The  fundamental  mode  Rayleigh  wave  is  exhib¬ 
ited  by  one  low  frequency,  low  velocity  zone.  The  first  higher  Rayleigh  mode  is 
exhibited  with  a  second  region  of  higher  frequency,  higher  velocity.  Figure  91  is 
the  phase  velocity  processing  of  the  35  meter  Betsy  data  set.  The  results  are 
very  complex.  This  complexity  continues  with  the  data  from  the  65  meter  offset 
(Rgure  92).  The  air  blast  cuts  through  the  data,  and  couples  with  the  higher 
mode  making  the  higher  mode  Rayleigh  waves  difficult  to  interpret.  The  funda¬ 
mental  mode  Rayleigh  wave  is  clear  at  the  5  meter  offset.  The  fundamental 
mode  has  a  limited  frequency  content  (to  25  Hz)  matching  the  results  of  the 
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Figure  87.  P-omega  stacking  of  5  meter  SWIG  data.  Nray  or  the 
number  of  velocities  computed  is  set  to  40. 
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Figure  88.  P-omega  stacking  o1 35  meter  SWIG  data.  Nray  or  the 
number  of  velocities  computed  is  set  to  40. 
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Phase  velocity  from  p-omega 
SWIG  data  from  5,35  &  65.0  meters 


Figure  89.  Interpreted  Love  wave  phase  velocity  dispersion  at  5, 
35  and  65  m.  5  meter  data  Is  represented  by  circles,  35  meter  by 
plus,  and  65  meter  by  triangles. 


PHASE  VELOCITY  (M/SEC)  *  10" 

0  00  1.00  2.00  3.00  4.00  5  00 


130 


.Sbei4  1st  DT  =0  001  SEC  NSAMP*1000 


’  - 1 - 1 - ! - 1 - 1  f 

.00  1  3.33  26.67  40.00  53.33  66.67  80. 


FREQUENCY  (HZ) 


Figure  90.  P-omega  array  processing  for  phase  velocity  of  5  meter 
Betsy  data.  Largest  amplitude  for  each  velocity  is  represented  by 
squares,  second  largest  by  circles,  third  largest  by  triangles,  and 
fourth  largest  by  plus  signs. 
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Figure  91 .  P-omega  array  processing  for  phase  velocity  of  35 
meter  Betsy  data.  Largest  amplitude  for  each  velocity  is 
represented  by  squares,  second  largest  by  circles,  third  largest  by 
triangles,  and  fourth  largest  by  plus  signs. 
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Figure  92.  P-omega  array  processing  for  phase  velocity  of  65 
meter  Betsy  data.  Largest  amplitude  for  each  velocity  is 
represented  by  squares,  second  largest  by  circles,  third  largest  by 
triangles,  and  fourth  largest  by  plus  signs. 
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group  velocity  processing.  The  energy  above  this  frequency  is  transferred  into 
the  first  higher  mode. 

The  interpretation  of  the  processed  results  was  more  difficult  than  for  the 
Love  wave  data.  The  fundamental  and  first  higher  modes  for  all  3  offsets  are 
given  in  Figure  93  and  appear  consistent  with  one  another.  As  noted  for  the 
Love  wave,  this  is  a  function  of  the  fact  that  the  phase  velocity  for  all  three  data 
sets  Imaged  the  same  Interval.  The  number  of  frequency  points  in  the  first 
higher  mode  data  set  was  so  small  as  to  limit  the  value  of  this  portion  of  the  data 
set  for  inversion  purposes. 
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Rayleigh  Wave  Fundamental  &  First  Higher  Mode 
Betsy  data  from  5.35  &  65  meters 
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Figure  93.  Interpreted  phase  velocity  dispersion  results  lor  the 
Betsy  data  set  at  5, 35  and  65  meter  offset.  Fundamental  mode  at 
5  meter  represented  by  plus,  at  35  meter  by  circle,  and  at  65  meter 
by  triangle.  First  higher  mode  phase  velocity  at  5  meter 
represented  by  triangle  hat,  at  35  meter  by  diamond  and  at  65  m 
byX. 


CHAPTER  4 

INVERSION  FOR  VELOCITY  STRUCTURE 


Much  effort  has  been  devoted  to  deriving  accurate  dispersion  curves  for 
the  surface  waves.  The  techniques  and  results  are  summarized  in  Chapter  3. 
These  dispersion  results  can  be  inverted  to  derive  a  velocity  and  Q  model.  As 
with  any  non-linear  inversion  technique,  a  starting  model  was  necessary.  In  the 
case  of  this  study  a  simple  refraction  analysis  of  body  wave  arrival  times  (both  P 
and  SH)  was  used  to  build  the  starting  model. 


1 .  Experimental  Technique 

The  direct  travel  path  time  to  any  receiver  can  be  computed  from  the 
source-receiver  spacing  and  the  surface  layer  velocity: 


T  =  — 
V,' 


where  T  is  the  total  travel  time,  xj  is  the  i'th  source  receiver  offset  and  V-j  is  the 
velocity  in  layer  1  or  the  surface  layer.  If  the  next  deeper  layer  has  a  velocity  V2 
that  is  greater  than  Vi ,  at  some  receiver  offset  the  first  arrival  or  fastest  travel 
time  will  be  from  the  deeper  layer.  The  travel  time  for  any  source  receiver  pair 
can  be  calculated  for  a  simple  layered  model  case.  Following  Telford  et  al. 
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(1976)  the  travel  times  for  a  2  layer  model  can  be  determined  from  the  three 
segments  of  the  travel  path  (Figure  94): 


(4.2) 


t  = 


OM 

V, 


MP  PR 


II  •  *u  I  *■  u-  sln0,  sine,  .  «  V, 

Using  the  relationship  -  --  '  =  or  sin  6q  =  — 

V,  V2  Vj 

metric  relationships,  equation  4,2  can  be  reduced  to: 


and  some  trigono- 


(4.3) 


2Z  cos  9 
V, 


This  equation  demonstrates  that  the  velocities  can  be  determined  from  the  slope 
of  the  first  arrivals  on  the  refraction  data.  The  intercept  time,  t,,  or  time  for  x=0, 
can  be  used  to  determine  the  layer  thickness  Z  from: 


(4.4) 


Z  = 


=  i  V,t, 


XOSB' 


This  analysis  can  be  extended  to  a  multiple  layer  model  (Figure  95)  assuming 
that  V-|  <  V2  <  V3.  The  refraction  path  OM'M"P"P'R’  is  fixed  according  to 
Snell's  law: 


sin  e,  sin 

V2 


sin  83 
V3 


(4.4) 


V, 


137 


where  02  is  the  critical  angle  for  the  second  layer  while  e-j  is  less  than  the  critical 
angle  for  the  upper  layer.  As  before,  the  travel  times  can  be  represented  as; 


(4.5) 


t  = 


OM'+R‘P*  M'M**+P’P" 

-  + - - - + 


V, 


V, 


M”P” 

V3 


As  with  the  two  layer  model,  using  some  trigonometric  relationships  equation  4.5 
can  be  rearranged  into: 


(4.6) 


cos  02  +  cos  0, 


This  equation  demonstrates  that  the  time-distance  curve  for  the  third  layer  has  a 
slope  inversely  proportional  to  the  velocity  of  the  third  layer  with  a  zero  intercept 
dependent  on  the  thickness  of  the  two  upper  layers.  This  conclusion  can  be 
generalized  to  any  number  of  layers  n,  such  that: 


(4.7)  t  =  +  2;^cos0„ 

where  0;  =  sin'’()^J.  TTiis  equation  can  be  used  for  any  number  of  layers 
providing  that  the  velocities  are  increasing  with  depth  and  that  each  layer  is  thick 
enough  to  permit  correct  analysis  from  the  time-depth  relationship.  The  individ¬ 
ual  velocities  and  the  subsequent  critical  angles  can  be  computed  from  the  slope 
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of  the  time-distance  curve.  The  layer  thicknesses  can  then  be  derived  from  the 
projected  intercepts  of  the  time-distance  curves. 

2.  Experimental  Results 

The  data  sets  for  this  study  were  acquired  from  5  to  94.5  meters  at  half 
meter  spacings.  All  data  sets  were  interpreted  using  the  Intercept  Time  Method 
(ITM)  as  documented  in  Telford  et  al  (1976).  The  basic  analysis  is  performed 
using  equation  4.7.  This  technique  assumes  non-dipping,  planar  layers  with  ve¬ 
locity  increasing  with  depth.  The  interpreted  first  break  picks  and  slopes  for  all 
six  data  sets  are  plotted  in  Figures  96  to  101.  The  data  is  plotted  on  an  ex¬ 
panded  scale,  both  in  time  and  distance.  The  data  is  summarized  in  tabular 
form  in  Tables  4  and  5.  Both  Betsy  and  SWIG  data  indicate  4  layers  with  similar 
thicknesses  for  the  corresponding  layers.  The  thicknesses  for  the  Betsy  data  set 
or  compressional  analysis  were  used  for  input  into  the  inversions.  The  primary 
reason  for  using  the  compressional  layer  thicknesses  was  that  the  inversion 
process  assumed  a  fixed  compressional  starting  model. 

B.  Inversion  of  Dispersion  Data 
1.  Experimental  Technique 

Detailed  explanations  of  the  inversion  algorithms  and  techniques  are  pre¬ 
sented  in  Herrmann  (1987).  Some  of  the  concepts  necessary  to  interpret  the  re¬ 
sults  of  this  thesis  are  presented  in  this  section.  Several  material  properties  in¬ 
fluence  surface  wave  dispersion,  including  shear  and  compressional  wave  ve¬ 
locity,  thicknesses  of  the  layers,  as  well  as  the  bulk  density  of  the  layers.  The  in¬ 
version  of  all  of  these  parameters  over  a  multilayer  model  is  a  nonlinear  process. 


Figure  97. 
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Figure  100.  Instrument  corrected  SWIG  data  at  35  m  with 
retraction  analysis  marked 
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Table  4 

BETSY  DATA  REFRACTION  ANALYSIS 


Layer 

Number 

1 

2 

3 

4 

P  Velocity 
(m/sec) 

408 

582 

666 

1143 

Intercept 

(msec) 

7.0 

12.5 

23.0 

57.5 

Crossover 

(meters) 

9.5 

57.0 

87.5 

NA 

Thickness 

(meters) 

1.6 

5.9 

10.6 

Halfspace 

Table  5 

SWIG  DATA  REFRACTION  ANALYSIS 


Layer 

Number 

1 

2 

3 

4 

SH  Velocity 
(m/sec) 

224 

337 

380 

482 

5.0 

19.0 

26.5 

65.7 

Crossover 

(meters) 

9.0 

23.0 

78.5 

NA 

Thickness 

(meters) 

2.1 

2.3 

10.6 

Halfspace 
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The  inversion  of  this  data  is  therefore  handled  in  an  iterative  fashion,  using  a 
damping  factor  to  limit  how  far  the  parameters  are  allowed  to  change  at  one 
time. 

One  additional  parameter  that  influences  the  velocity  dispersion  is  the  at¬ 
tenuation  or  Q  factor  of  any  layer  or  material  relative  to  compressional  or  shear 
wave  propagation.  The  attenuation  of  the  surface  wave  due  to  Q  must  be  taken 
into  account  in  order  to  keep  the  results  causal.  The  first  step  in  the  inversion 
processing  was  to  do  a  non-causal  inversion,  or  an  inversion  that  neglected  the 
effects  of  Q.  The  primary  reason  for  doing  a  non-causal  inversion  was  that  the 
attenuation  appeared  to  be  poorly  measured,  and  therefore  it  seemed  appropri¬ 
ate  to  obtain  a  starting  result  that  was  not  contaminated  by  the  uncertainty  in 
measuring  attenuation  due  to  Q.  Once  a  non-causal  result  was  obtained,  a 
causal  inversion,  or  inversion  in  which  Q  influenced  dispersion  velocity  could  be 
performed. 

Because  there  are  so  many  parameters  influencing  the  dispersion,  certain 
assumptions  about  the  material  properties  were  made  to  constrain  the  final  re¬ 
sults.  One  assumption  that  was  input  into  the  inversion  was  that  the  shear  Q  or 
Op  was  related  to  the  compressional  Q  or  by  a  constant,  in  this  case  2.25. 

For  this  study  it  was  also  assumed  that  the  number  of  layers,  layer  thicknesses 
and  compressional  velocities  were  all  known  from  the  refraction  analysis.  This 
leaves  velocity  dispersion  and  surface  wave  attenuation  as  the  observations  to 
invert  to  change  the  shear  velocity  and  shear  Q  for  the  various  layers.  These 

data  can  be  inverted  simultaneously  with  non-zero  partial  derivatives  for  attenu¬ 
ation  versus  p  and  dispersion  velocity  versus  Qjj.  This  type  of  inversion  is 


149 


known  as  a  coupled  inversion.  A  decoupled  inversion  is  where  either  Qp  or  |3 

are  assumed  constant  and  the  other  parameter  is  determined  by  inversion. 

In  the  processing  of  this  data,  decoupled  inversions  were  used  to  obtain 
estimates  of  the  parameter  values  before  performing  a  final  coupled  inversion. 
Equation  1.10  demonstrates  that  the  phase  velocity  and  group  velocity  of  a  sur¬ 
face  wave  in  a  plane  layered  medium  are  directly  related  to  one  another.  Since 
one  velocity  can  be  computed  from  the  other,  the  two  data  sets  should  be 
equivalent  for  inversion  purposes.  The  inversions  that  were  performed  included 
using  the  group  and  phase  velocities  separately  as  well  as  together  to  solve  for 
the  Q|3  and  (3. 


2.  Experimental  Results 

The  inversion  software  as  implemented  by  Herrmann  (1987),  includes  the 
capability  to  invert  for  Q  structure  as  well  as  velocity  structures.  Inputs  can  be 
any  combination  of  Love  group  velocity,  phase  velocity,  and  gamma  or  attenu¬ 
ation  coefficients  as  well  as  the  corresponding  Rayleigh  characteristics.  The 
Rayleigh  and  Love  wave  information  were  inverted  as  separate  problems  in  this 
study.  The  two  data  sets  may  be  imaging  different  information,  such  as  SV  and 
SH.  Additionally,  the  Rayleigh  information  was  more  complex  and  warranted 
separate  treatment. 

The  inversion  techniques  can,  in  principle,  invert  multiple  modes  of  veloc¬ 
ity  information.  One  limitation  is  that  these  different  modes  must  encompass  the 
same  frequency  band.  The  energy  in  the  frequency  band  below  25  Hz  was  pri¬ 
marily  in  the  fundamental  mode  while  energy  from  the  frequencies  above  25  Hz 
were  in  the  first  higher  mode.  This  result  is  similar  to  the  Simila  (1982)  and 
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Bogaards  (1989)  study.  Since  these  two  modes  did  not  overlap  in  frequency, 
they  could  not  be  inverted  together.  Additionally,  the  algorithms  for  different 
modes  would  not  allow  for  the  inversion  of  a  higher  mode  by  itself.  Conse¬ 
quently,  inversion  of  the  Rayleigh  wave  data  was  limited  to  the  fundamental 
mode  information. 

One  Rayleigh  wave  inversion  was  performed  for  this  thesis,  utilizing  the 
maximum  amount  of  available  data.  The  phase  velocity  from  the  p-omega  array 
technique  was  combined  with  the  group  velocity  and  attenuation  coefficients 
from  the  Barker  array  dispersion  analysis.  The  starting  model  included  com- 
pressional  and  shear  velocities  from  the  refraction  analysis  and  Q  values  from 
the  Grant  (1988)  model  (Table  3).  The  inversion  was  run  as  a  causal  inversion, 
where  Q  and  velocity  structure  both  affect  the  surface  wave  dispersion.  The  re¬ 
sults  of  the  velocity  inversion  are  displayed  in  Figure  102  with  velocity  versus 
depth  along  with  the  resolution  kernels.  The  normalized  amplitude  of  the  resolu¬ 
tion  kernels  for  each  layer  are  presented  with  depth  to  the  center  of  the  layer 
across  the  top  versus  the  depth  for  each  value  of  the  kernels.  The  results  dem¬ 
onstrate  that  the  kernels  have  maximum  amplitudes  for  the  appropriate  layers. 
The  unnormalized  kernels  indicate  that  the  amplitudes  for  these  kernels  de¬ 
crease  with  depth.  The  corresponding  final  decoupled  inversion  for  Q  with 
resolution  kernels  is  displayed  in  Figure  103.  The  third  model  in  Figures  102 
and  103  is  the  coupied,  composite  inversion  of  the  velocity  information  with  the 
Q  or  attenuation  information.  In  the  composite  inversions  the  velocity  and  at¬ 
tenuation  information  are  normalized  since  the  amplitudes  may  be  of  entirely  dif¬ 
ferent  orders  of  magnitude.  This  allows  the  inversion  to  not  be  overly  weighted 
or  biased  by  either  the  Q  or  the  velocity  information.  Figure  1 04  shows  the  fit 
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Rayleigh  wave  65m  shear  velocity  inversion 
SIIF  AR  Vtl.OCITY  (  M/SF.Cl  RESOLVING  KERNELS 


Rgure  102.  65m  Betsy  inversion  for  shear  velocity  structure. 
Diamonds  represent  the  start  or  input  model;  triangles  are  the 
decoupled  causal  inversion;  squares  are  the  final  coupled,  causal 
inversion  results.  Horizontal  bars  are  standard  deviation. 
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Sun  Model  with  Qb  inveriion 
Q  BETA  INV  RESOLVING  KERNELS 

0.06  0.12  0.18  0.24  NORMALIZED 


txqbsien.mod 
dtq  b  1  .  mod 
□  qb  3 . mod 


Rgure  103.  65m  Betsy  inversion  for  shear  Q  inverse  structure. 
X,s  represent  the  start  or  input  model;  triangles  are  the  decoupled 
causal  inversion;  squares  are  the  final  coupled,  causal  inversion 
results.  Horizontal  bars  are  standard  deviation. 
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Gamma  Output  from  coupled  causal  inversion 
Using  Betsy  65m  data 

0  090  - - 


0.07  5  -I 


I 

! 


0.05  9 


0.044  J 


0.029 


0.014 


i  ' 

/ 


¥ 


IRF  ^ 


-0.002 


^  disp.gamSf 


15 


30  45  60 

FREQUENCY  (HZ) 


RAY'LEIGH  WAVE  | 
75  90 


Figure  104.  Comparison  of  observed  to  predicted  gamma  for  65m 
Betsy  inversion.  Triangles  represent  observed  data  with  vertical 
bars  representing  one  standard  deviation.  The  solid  curve  is  the 
predicted  results. 
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between  the  observed  attenuation  coefficients  and  the  predicted  attenuation  co¬ 
efficients  from  the  final  inverted  velocity  and  Q  models.  Figure  1 05  shows  the  fit 
between  the  observed  and  predicted  group  and  phase  velocities.  The  fit  is  much 
better  for  the  group  velocities  than  the  phase  velocities  since  the  group  velocity 
has  many  more  points  in  the  data  set  and  the  inversion  is  trying  to  minimize  the 
least  squares  error  between  the  predicted  and  observed  data  sets.  The 
inconsistency  in  the  results  between  the  group  and  phase  velocities  may  reflect 
the  fact  that  the  phase  velocities  and  group  velocities  were  imaging  different 
lateral  extents  of  the  velocity  structures. 

The  inversion  process  was  applied  to  the  Love  waves  utilizing  three  dif¬ 
ferent  portions  of  the  data  set.  The  phase  velocity  and  gamma  or  attenuation 
coefficient  data  were  inverted  first.  The  group  velocity  and  gamma  were  inverted 
second.  The  final  inversion  included  the  group  velocity,  the  phase  velocity  and 
the  gamma  data.  Each  of  these  inversions  followed  the  same  three  stages.  The 
first  was  to  invert  for  a  non-causal  velocity  field,  or  an  inversion  for  shear  velocity 
structure  that  ignored  the  effects  of  Q.  The  second  inversion  included  a  decou¬ 
pled,  causal  inversion  for  shear  velocity  structure.  In  this  instance  the  Q  struc¬ 
ture  was  held  constant  while  allowing  the  shear  velocity  to  vary.  The  third  inver¬ 
sion  was  a  coupled,  causal  inversion  which  allowed  both  the  shear  velocity  struc¬ 
ture  to  vary  based  on  the  dispersion  inputs,  and  the  Q  structure  to  vary  based  on 
the  gamma  inputs. 

Figure  106  summarizes  the  results  of  the  Love  wave  phase  velocity  in¬ 
version  for  shear  velocity.  The  shear  velocities  in  the  shallowest  two  layers  cal¬ 
culated  from  the  refraction  analysis  appear  to  be  too  high.  There  are  several 
possible  explanations  for  this  result.  One  is  that  there  could  be  some 
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experimental  error  since  there  are  only  a  few  time  depth  pairs  used  in  the  shal¬ 
lowest  layers.  A  more  reasonable  explanation  would  be  that  the  dispersion  re¬ 
sults  are  averaged  over  a  longer  spatial  interval,  on  the  order  of  30  meters  for 
the  phase  velocity  and  at  least  65  meters  for  the  group  velocity.  The  refraction 
analysis  for  the  shallowest  layer  covered  only  9  meters  laterally.  The  two  deep¬ 
est  layers  appear  to  be  fairly  stable  over  all  four  sets  of  inversions  (three  Love 
and  one  Rayleigh  wave).  The  maximum  of  the  resolution  kernels  for  both  of  the 
deepest  layers  is  about  10*2  which  indicates  that  the  longest  wavelengths  are 
not  contributing  much  to  the  model.  Dispersion  data  below  the  6  Hz  cutoff  in  this 
data  set  would  help  resolve  these  layers. 

Figure  107  illustrates  the  Q  structure  from  the  phase  velocity  inversion. 
The  results  show  significant  scatter  for  the  different  inversions,  and  probably  can 
only  yield  a  rough  estimate  of  Q.  The  primary  input  into  this  inversion  is  a  com¬ 
parison  of  the  spectral  amplitudes  with  range  for  the  Love  wave  portion  of  the 
data  set.  The  assumption  implicit  in  this  technique  is  that  spreading  loss  and  Q 
are  the  only  parameters  influencing  relative  amplitude  of  the  Love  waves.  Some 
spectral  notching  due  to  geology  was  noted  during  the  section  on  MFT.  This 
could  very  easily  contaminate  the  gamma  data  and  bias  the  Q  inversion  results. 

The  predicted  phase  velocity  dispersion  from  the  final  inverted  velocity 
model  is  given  in  Figure  108.  The  visual  fit  is  very  good,  in  spite  of  the  fact  that 
the  input  data  have  fairly  large  error  bars,  especially  at  the  low  frequencies.  The 
fit  of  the  predicted  to  observed  gamma  data  is  shown  in  Figure  109.  The  fit  is 
fair  with  most  of  the  divergence  at  the  low  end. 

A  similar  con  parison  of  the  group  velocity  results  is  shown  in  Figures  110 
to  1 14.  The  anomalous  thing  about  these  results  is  that  after  the  Q  inversion. 
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Velocity  Output  from  coupled  causal  inversion 
Betsy  65m  data  used  as  input 
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Figure  105.  Comparison  of  observed  to  predicted  group  and  phase 
velocities  for  65m  Betsy  inversion(Rayleigh  surface  wave).  The 
RPF  triangles  represent  observed  phase  velocity  data  and  the  RGF 
triangles  represent  the  group  velocity.  The  solid  curve  in  both 
cases  is  the  predicted  result. 
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l.ovc  wave  65fn  C  inversion 

veil  non^causal.  vc)2.  decoupled:  vel3  coupled 
SHEAR  VELOCITY  (  M/SFC)  RESOLVING  KERNELS 


Figure  106.  65  meter  SWIG  data  phase  velocity/gamma  inversion 
tor  shear  velocity  structure.  Plus  is  input  start  model,  squares  is 
non-causal  inversion,  circles  is  decoupled,  causal  inversion  with 
fixed  Q  model,  triangles  is  final  model  from  coupled  causal 
inversion. 
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Figure  107.  65  meter  SWIG  data  phase  velocity/gamma  inversion 
tor  Q  structure.  Plus  is  input  start  model,  squares  is  final  model 
from  coupled  causal  inversion. 
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Figure  108.  65  meter  SWIG  data  observed  versus  predicted  phase 
velocity  after  final  inversion.  Triangles  are  observed  data  with 
vertical  bars  being  standard  deviation  from  array  processing 
analysis.  Solid  line  is  predicted  data  set. 


GAMMA  (1/M) 


160 


0.060  1 

0.045  ■ 

0.030  - 

I 

j 

0.014  -j 

I 

:  i  i 

LF-ii^  disp.gamSf  LOVE  wave  j 

-0.001  -i - - 1 - - - - - 

0  15  30  45  60  75  90 


FREQUENCY  (HZ) 


Rgure  109.  65  meter  SWIG  data  observed  versus  predicted 
attenuation  coetlicients  after  final  inversion.  Triangles  are 
observed  data  with  vertical  bars  being  standard  deviation  from 
array  processing  analysis.  Solid  line  is  predicted  data  set. 


I 

ff 


DKPTH  (Ml 


161 


Love  wave  65in  I'  inversion 

veil:  non-causti;  vel2.  decoupled:  vel3:  coupled 
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Figure  110.  65  meter  SWIG  data  group  velocity/gamma  inversion 
for  shear  velocity  structure.  Plus  is  input  start  model,  squares  is 
non-causal  inversion,  circles  is  decoupled,  causal  inversion  with 
fixed  O  model,  triangles  is  final  model  from  coupled  causal 
inversion. 
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Figure  111.  65  meter  SWIG  data  group  velocity/gamma  inversion 
for  Q  structure.  Plus  is  input  start  model,  squares  is  final  model 
from  coupled  causal  inversion. 
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Figure  112.  65  meter  SWIG  data  observed  versus  predicted  group 
velocity  after  second  inversion.  Triangles  are  observed  data  with 
vertical  bars  being  standard  deviation  from  array  processing 
analysis.  Solid  line  is  predicted  data  set. 
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Love  wave  65m  U  inversion 
vel3:  coupled 
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Figure  113.  65  meter  SWIG  data  observed  versus  predicted  group 
velocity  after  third  and  final  inversion.  Triangles  are  observed  data 
with  vertical  bars  being  standard  deviation  from  array  processing 
analysis.  Solid  line  is  predicted  data  set. 
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Figure  114.  65  meter  SWIG  data  observed  versus  predicted 
attenuation  coefficients  after  final  inversion.  Triangles  are 
observed  data  with  vertical  bars  being  standard  deviation  from 
array  processing  analysis.  Solid  line  is  predicted  data  set. 
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the  fit  to  the  group  velocity  curve  is  degraded.  The  excellent  comparison  of  the 
predicted  and  observed  group  velocity  after  inversion  with  the  Q  held  constant  is 
shown  in  Figure  112.  Allowing  the  Q  to  vary  in  the  coupled  causal  inversion 
again  resulted  in  poorer  fits  to  both  group  velocityfFigure  113)  and 
gamma(Figure  114).  The  degradation  of  results  implies  some  inconsistency  be¬ 
tween  the  different  portions  of  data.  The  most  probable  difficulty  is  that  the  at¬ 
tenuation  coefficients  are  not  strictly  a  function  of  Q. 

The  final  set  of  inversions  incorporates  all  three  sets  of  data,  group  veloc¬ 
ity,  phase  velocity  and  gamma.  The  inverted  results  are  included  in  Figures  115 
to  118.  This  analysis  includes  the  maximum  amount  of  input  data  and  should 
therefore  provide  the  best  results.  Figure  117  demonstrates  that  the  final  model 

is  consistent  with  both  the  input  group  and  phase  velocity,  although  (Figure  118) 
the  final  model  did  not  produce  a  physical  (greater  than  zero)  Qp  structure  that 

was  consistent  with  the  measured  attenuation  coefficients.  The  variability  of  in¬ 
verted  Q|5  and  discrepancies  with  predicted  attenuation  coefficients  indicates 

that  this  technique  will  only  yield  a  rough  estimate  of  Q.  Additionally,  it  indicates 
that  the  attenuation  with  frequency  is  influenced  by  more  factors  than  the  two 
assumed  in  this  study,  simple  spreading  loss  and  loss  due  to  formation  Q. 

C.  Forv/ard  Modeling 

Previous  dispersion  analysis  including  Simila  (1982)  have  used  forward 
modeling  techniques  to  match  a  velocity  structure  to  measured  velocity  disper¬ 
sion.  In  general,  an  attempt  was  made  to  match  observed  to  predicted  disper¬ 
sion.  In  this  study,  an  attempt  will  be  made  to  match  the  entire  synthetic  wave¬ 
form  to  the  measured  waveforms.  Because  there  can  be  more  than  one  model 
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Figure  115.  65  meter  SWIG  data  combined  group  and  phase 
velocity/gamma  inversion  for  shear  velocity  structure.  Plus  is  input 
start  model,  squares  is  non-causal  inversion,  circles  is  decoupled, 
causal  inversion  with  fixed  Q  model,  triangles  is  final  model  from 
coupled  causal  inversion. 
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Figure  117.  65  meter  SWIG  data  observed  versus  predicted  group 
and  phase  velocity  after  final  inversion.  Triangles  are  observed 
data  with  vertical  bars  being  standard  deviation  from  array 
processing  analysis.  Solid  line  is  predicted  data  set. 
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Figure  118.  65  meter  SWIG  data  observed  versus  predicted 
attenuation  coefficients  after  final  inversion  of  group  velocity,  phase 
velocity  and  gamma.  Triangles  are  observed  data  with  vertical 
bars  being  standard  deviation  from  array  processing  analysis. 

Solid  line  is  predicted  data  set. 
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Table  6 

INPUT  PARAMETERS  FOR  SYNTHETIC  SURFACE 
WAVE  COMPUTATIONS 


Parameter 

Layer  1 

Layer  2 

Layer  3 

Layer  4 

Layer 

Thickness 

1.6 

5.9 

10.6 

Halfspace 

Qa 

9.9 

22.5 

22.5 

22.5 

Qp 

4.4 

10 

10 

10 

a 

408 

582 

666 

1143 

P 

146.7 

250.7 

424.1 

619.8 

e 

1.8 

1.9 

1.9 

2.0 
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that  matches  a  given  dispersion  data  set,  a  good  visual  match  of  synthetic  data 
would  add  confidence  to  the  calculated  velocity  model. 

A  synthetic  Love  wave  was  calculated  for  the  model  in  Table  6,  which  was 
the  velocity  results  from  the  final  Love  wave  inversion  combined  with  Q  values 
from  the  Grant  (1988)  study  which  appeared  more  physically  reasonable  than 
the  ones  computed  for  this  study.  The  algorithms  for  computing  synthetic  sur¬ 
face  waves  are  fully  documented  in  Herrmann  (1985).  The  results  are  compared 
to  the  Love  wave  stripped  from  the  data  set  using  phase  matched  filtering  in 
Chapter  3  and  displayed  in  Figure  77.  The  first  difficulty  encountered  was  that 
the  synthetic  data  sets  did  not  match  the  acquired  data  set  in  frequency  content. 
The  spectral  content  of  the  synthetics  was  flat  below  the  center  frequency  (about 
10  Hz)  and  dropped  off  sharply  above  the  center  frequency.  A  triangular  source 
function  with  a  40  msec  width  was  finally  used  as  a  best  match.  The  synthetics 
appeared  to  have  the  opposite  polarity  to  the  real  data,  so  all  results  have  been 
flipped  for  a  better  match.  The  unfiltered  output  is  presented  in  Figure  119  on  a 
scale  comparable  to  the  PMF  output  in  Figure  77.  The  synthetic  Love  wave 
does  not  match  as  well  as  would  be  desired,  with  the  largest  difference  being 
apparent  in  the  number  of  lobes  and  apparent  frequency  content.  The  low  end 
of  the  spectrum  predominates  in  this  data  set  relative  to  the  sharp  rolloff  in  the 
high  frequency  end.  A  high  pass  filter  at  10  Hz  was  applied  to  the  data  set  and 
is  displayed  as  Figure  1 20.  A  better  fit  was  obtained  by  filtering  at  20  Hz  and  is 
presented  in  Figure  121.  The  PMF  results  have  a  high  frequency  ripple  at  the 
end  corresponding  to  slower  velocities.  It  appears  that  these  high  frequency  re¬ 
sults  from  the  PMF  may  not  be  real,  in  that  the  input  constraints  to  the  PMF  were 
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limited  to  50  Hz.  Data  outside  this  frequency  limit  would  not  have  been  elimi¬ 
nated  with  a  boxcar  type  filter  due  to  the  problem  with  Gibb’s  phenomena. 

Several  attempts  were  made  to  boost  the  frequency  content  of  the  syn¬ 
thetic  with  little  success.  When  the  triangular  source  function  was  cut  back  from 
40  msec,  the  data  exhibited  a  high  frequency  ripple  on  the  early  arrival  part  of 
the  wave  train.  An  example  of  this  is  Figure  122,  where  the  data  has  had  a  tri¬ 
angular  function  of  8  msec  with  a  low  cut  filter  of  20  Hz  applied.  It  seems  sur¬ 
prising  that  the  high  frequency  portion  of  this  data  set  arrives  earlier  or  faster  in 
the  data  set,  since  the  dispersion  analysis  demonstrates  that  the  highest  fre¬ 
quencies  should  come  in  latest. 

Several  Rayleigh  wave  synthetics  were  also  generated.  Because  there 
were  two  modes,  the  synthetics  were  computed  three  ways.  The  first  synthetic 
in  Figure  123  contains  the  fundamental  mode  only.  Figure  124  contains  the  first 
higher  mode  and  Figure  125  is  the  results  of  a  synthetic  containing  both  the  fun¬ 
damental  as  well  as  the  first  higher  mode  Rayleigh  wave.  Some  of  the  more  im¬ 
portant  features  of  the  real  data  are  not  replicated  in  the  synthetic.  The  two 
Rayleigh  modes  in  the  data  are  composed  of  non-overlapping  frequencies.  It  did 
not  appear  possible  to  generate  such  a  bimodal  synthetic,  matching  the  appro¬ 
priate  spectral  contents. 
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Figure  120.  Love  wave  synthetic  with  low  cut  filter  at  10  Hz. 
Triangular  source  width  40  msec. 
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Figure  121 .  Love  wave  synthetic  with  low  cut  filter  at  20  Hz, 
Triangular  source  width  40  msec. 
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Figure  122.  Love  wave  synthetic  with  low  cut  filter  at  20  Hz, 
Triangular  source  width  8  msec. 


Figure  1 23.  Rayleigh  wave  synthetic  of  fundamental  mode  only 
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Figure  124.  Rayleigh  wave  synthetic  of  first  higher  mode  only 
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Figure  1 25.  Rayleigh  wave  synthetic  of  combined  fundamental  and 
first  higher  mode 


CHAPTER  5 
CONCLUSIONS 


The  objective  of  this  study  was  the  determination  of  the  shear  velocity 
structure  in  a  shallow  alluvial  environment.  This  thesis  utilized  several  tech¬ 
niques  that  had  been  used  individually  in  previous  work  applied  to  different  scale 
problems.  The  project  was  composed  of  three  major  segments.  The  first  seg¬ 
ment  was  the  acquisition  of  the  field  data.  Second,  the  field  data  were  proc¬ 
essed  to  determine  dispersion  curves.  The  third  segment  was  calculating  a  ve¬ 
locity  model  by  inverting  the  dispersion  information.  This  final  section  will  try  to 
bring  together  the  different  techniques  and  propose  how  they  can  b:  used  in  fu¬ 
ture  work. 

A.  Acquisition  Techniques 

The  field  data  were  acquired  with  an  array  of  60  receivers  spaced  at  0.5 
meter  increments.  The  sources  were  set  at  3  ranges  from  the  end  of  the  re¬ 
ceiver  array,  at  5.0  meters,  at  35.0  meters  and  at  65.0  meters.  Two  sources 
were  used  in  this  study,  the  Betsy  and  the  SWIG.  The  Betsy  source  generated  P 
wave  data  as  well  as  a  Rayleigh  surface  waves.  The  SWIG  source  provided  SH 
energy  and  therefore  Love  waves.  The  results  for  the  SWIG  were  always  the 
most  easily  interpreted.  The  Betsy  data  set  was  more  difficult  to  interpret  for 
three  reasons:  first,  the  Betsy  data  contained  two  modes  of  Rayleigh  waves  and 
the  SWIG  data  set  only  contained  one  surface  wave  mode;  second,  the  Betsy 
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data  set  contained  an  airblast  arrival  which  actively  coupled  with  the  first  higher 
mode  Rayleigh:  and  finally,  the  time  separation  between  the  body  waves  and  the 
surface  waves  was  greater  for  the  SWIG  data  set. 

The  airblast  biased  the  first  higher  mode,  making  it  useless  for  inversion 
purposes.  The  fundamental  Rayleigh  mode  had  a  bandwidth  of  approximately  6 
to  26  Hz.  Almost  all  of  the  energy  above  26  Hz  appeared  in  the  first  higher 
mode  Rayleigh.  In  theory,  having  two  modes  of  a  surface  wave  should  help 
constrain  the  inversion,  yielding  a  final  answer  that  is  less  ambiguous.  This  extra 
data  could  not  be  used  in  this  study  since  the  inversion  technique  used  would 
only  incorporate  different  data  types  over  the  same  frequency  interval.  An  addi¬ 
tional  limitation  was  that  the  higher  mode  data  could  not  be  used  by  itself.  The 
Simila(1982)  data  acquired  with  a  different  source,  high  explosives,  at  McCor¬ 
mick  ranch  showed  a  similar  frequency  content  for  the  two  Rayleigh  wave 
modes.  This  demonstrates  that  a  different  source  would  not  alleviate  these 
problems  in  this  alluvial  environment. 

The  SWIG  data  remained  simple  and  easily  interpretable  through  the  en¬ 
tire  range  of  5  to  94.5  meters  used  in  this  study.  The  data  contained  a  large 
amplitude  direct  and  refracted  first  arrival.  These  events  were  followed  by  a  vis¬ 
ibly  dispersive,  single  mode  of  Love  wave.  The  time  separation  between  these 
two  portions  of  the  data  set  was  larger  for  the  SWIG  data  than  for  the  Betsy  data, 
which  made  the  dispersion  analysis  easier  to  perform. 

B.  Dispersion  Analysis 

Dispersion  analysis  of  the  field  data  yielded  group  velocity,  phase  velocity 
and  attenuation  information.  The  group  velocity  effectively  images  or  illuminates 


183 


the  average  structure  of  the  entire  interval  between  the  source  and  the  receiver. 
The  phase  velocity  is  a  function  of  the  average  structure  between  the  receivers 
used  in  the  determination  of  the  phase  velocity.  The  attenuation  coefficients 
were  based  on  evaluation  of  changes  in  spectral  amplitude  with  range  after  cor¬ 
rection  for  spreading  loss.  The  assumption  used  in  this  study  was  that  any  am¬ 
plitude  losses  were  due  to  the  attenuation  or  the  average  Q  of  the  formations  be¬ 
tween  the  source  and  the  receiver.  Inconsistency  in  Q  estimates  indicates  that 
the  simple  spreading  corrections  may  be  inappropriate. 

The  most  basic  dispersion  analysis  technique  used  was  the  multiple  filter 
analysis  (MFA)  which  determines  the  group  velocity  from  a  single  source  re¬ 
ceiver  pair.  The  MFA  computes  the  peaks  of  the  envelope  function  after  multiple 
narrow  band  Gaussian  filters  have  been  applied  to  the  single  trace.  These 
peaks  are  then  related  to  velocity  and  the  center  frequency  of  the  filters  to  derive 
a  dispersion  curve.  This  technique  is  computationally  less  intense  than  the  other 
techniques  used  in  this  study.  The  composite  display  output  from  this  analysis 
includes  a  plot  of  the  single  trace  versus  time  as  well  as  versus  velocity  (see 
Figures  25  to  45  and  Appendix  A).  This  information  along  with  the  spectral  am¬ 
plitude  plots  greatly  aids  the  interpreter  in  determining  what  is  going  on  in  the 
data  set.  This  analysis  can  enable  the  interpreter  to  determine  where  the  sur¬ 
face  waves  are,  and  what  modes  are  present.  One  very  important  use  of  the 
MFA  in  this  study  was  to  observe  changes  in  group  velocity  with  range.  These 
observations  can  be  useful  in  determining  at  what  range  the  surface  waves  have 
stabilized,  and  additionally,  to  determine  if  there  are  some  lateral  velocity  or 
structural  changes  over  the  range  of  the  data  set  such  as  was  shown  in  Figure 
46. 
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The  array  or  Barker  technique  essentially  averages  the  group  velocities 
from  multiple  source-receiver  pairs.  Additionally,  this  technique  yielded  the  at¬ 
tenuation  values  subsequently  used  in  the  inversion  process.  One  useful  feature 
of  this  array  technique  is  that  it  is  not  sensitive  to  the  spatial  sampling  interval.  In 
this  study,  one  bad  trace  in  the  data  set  could  be  thrown  out  prior  to  the  process¬ 
ing.  This  technique  yielded  very  consistent  results  for  both  the  Love  and 
Rayleigh  wave  data  sets. 

The  phase  matched  filter  (PMF)  technique  was  the  only  other  single  sta¬ 
tion  technique  utilized  in  this  study.  It  can  be  used  to  fine  tune  a  single  station 
answer  derived  from  some  other  technique,  which  for  this  study  was  the  MFA. 
The  technique  yielded  a  visually  better  answer  in  most  cases,  but  was  very  sen¬ 
sitive  to  interfering  events  closely  spaced  in  time.  One  additional  difficulty  with 
this  procedure  was  that  it  was  interactive  and  very  sensitive  to  user  inputs.  The 
final  group  velocities  were  very  similar  to  the  results  from  the  array  or  Barker 
technique.  Since  the  Barker  technique  was  designed  for  batch  processing  and  is 
relatively  insensitive  to  user  inputs,  the  PMF  technique  would  be  a  secondary 
choice  for  dispersion  analysis  if  array  data  was  available.  One  use  for  the  PMF 
would  be  as  a  filter.  Figure  77  shows  the  result  of  stripping  out  a  data  set  con¬ 
taining  only  the  fundamental  mode  Love  wave.  Alternatively,  this  technique 
could  be  used  to  remove  the  surface  wave  for  more  detailed  analysis  of  the  re¬ 
maining  body  waves. 

The  phase  velocity  was  obtained  from  the  p-omega  technique  which  is 
based  on  a  concept  provided  by  McMechan  and  Yedlin  (1981).  A  i-p  transform 
followed  by  a  1  -d  FFT  yielded  data  in  the  form  of  slowness  (p)  versus  frequency 
(omega).  Since  slowness  is  the  inverse  of  velocity,  the  results  can  be  plotted  as 
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phase  velocity  versus  frequency  for  fairly  straightforv/ard  interpretation.  The 
composite  results  for  all  three  SWIG  or  Love  wave  data  sets  are  shown  in  Figure 
89.  The  phase  velocities  for  all  3  sets  are  extremely  consistent.  The  Rayleigh 
data  sets  show  the  same  consistency  of  phase  velocity  for  different  ranges 
(Figure  93).  This  result  is  in  marked  contrast  to  the  group  velocity  results  where 
velocity  increases  substantially  with  range  (see  Figures  52,  66  and  67  for  group 
velocity  results). 

Equation  1.10  demonstrates  that  in  a  plane  layered  medium  it  is  possible 
to  calculate  the  phase  velocity  from  the  group  velocity  or  the  group  velocity  from 
the  phase  velocity.  The  phase  velocity  remained  consistent  with  range.  All  three 
SWIG  data  sets  had  almost  identical  phase  velocities,  however  the  group  veloc¬ 
ity  varied  with  range,  which  is  inconsistent  with  equation  1.10.  The  geometry  of 
the  experiment  (as  demonstrated  in  Figure  9)  indicates  that  the  group  velocity 
would  be  a  function  of  the  material  properties  of  a  larger  region  as  the  source 
offset  increases.  The  phase  velocity  is  a  function  of  the  material  properties  be¬ 
tween  the  receivers  in  the  array.  Since  the  receivers  did  not  move  in  this  study, 
the  phase  velocity  imaged  the  same  material  for  all  3  source  positions.  The  ap¬ 
parent  lateral  variation  in  material  properties  violates  the  plane  layered  assump¬ 
tion  for  the  inversion  techniques  that  were  used  in  this  study.  The  interval  inter¬ 
rogated  by  the  surface  waves  allowed  for  the  determination  of  the  average  ma¬ 
terial  properties. 


C.  Inversion  Technique 

The  inversion  portion  of  the  study  consisted  of  three  discrete  phases.  The 
first  phase  involved  normal  refraction  analysis  of  both  P  and  SH  data.  Next,  the 
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refraction  model  was  used  as  the  start  model  for  a  non-linear  inversion  of  the 
dispersion  analysis.  Finally,  synthetics  were  generated  from  the  inverted  model 
structure  for  comparison  to  the  original  field  data. 

The  refraction  analysis  was  performed  on  a  single  linear  array  assuming 
horizontal  layers.  The  technique  used  is  documented  in  Telford  et  al.  (1976)  and 
the  results  are  presented  in  Tables  4  and  5.  The  SWIG  data  set  provided  an  un¬ 
contaminated  SH  wavefield  for  direct  and  refracted  arrivals.  The  SH  results  ap¬ 
peared  consistent  with  the  compressional  refraction  analysis  except  for  some  dif¬ 
ferences  in  the  layer  thicknesses.  The  refraction  results  were  used  as  a  starting 
model  for  the  surface  wave  inversions  that  followed. 

The  Love  wave  inversions  yielded  shear  velocities  that  were  similar  to  the 
refraction  analysis  results  for  the  two  deepest  layers.  The  inversions  indicated 
the  shallowest  two  layers  were  slower  than  the  refraction  survey  results.  The 
shallow  layer  discrepancies  are  probably  a  result  of  the  fact  that  the  surface 
waves  imaged  a  much  longer  interval  than  the  refraction  survey  for  the  first  two 
layers.  The  resolution  kernels  were  very  small  for  the  deepest  two  layers,  which 
indicates  that  the  wavelengths  were  not  long  enough  to  image  at  these  depths. 
It  seems  obvious  that  data  at  a  lower  frequency  than  the  6  Hz  used  in  this  study 
would  have  been  beneficial.  The  other  significant  result  was  that  the  Q  inversion 
appeared  fairly  unstable,  frequently  trying  to  yield  negative  Q  values,  even  when 
the  models  were  heavily  damped.  This  instability  is  less  due  to  technique  than 
the  fact  that  the  input  attenuation  coefficients  are  assumed  to  only  be  a  function 
of  Q  and  spreading  loss.  There  appear  to  be  more  factors,  including  scattering, 
that  are  adding  to  the  complexity  of  the  attenuation  with  frequency. 
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The  final  step  of  this  project  was  to  generate  synthetic  surface  waves  to 
match  the  original  field  data  from  the  final  model.  This  synthetic  was  at  least 
partially  successful  for  the  Love  waves.  The  final  synthetic  wavefield  (Figure 
121)  has  a  frequency  content  and  character  similar  to  the  PMF  extracted  Love 
waves  in  Figure  77.  The  complexity  of  the  Rayleigh  wave  data  set  appeared  to 
be  beyond  the  capability  of  the  techniques  used  in  this  project.  The  biggest  diffi¬ 
culty  was  in  duplicating  the  behavior  of  the  splitting  of  the  energy  spectrum  be¬ 
tween  the  fundamental  and  first  higher  mode. 

The  final  velocity  model  determined  for  the  McCormick  Ranch  area  is  the 
model  that  was  input  into  the  synthetics  and  documented  in  Table  6.  The  layer 
thickness  for  the  shallowest  layer  is  given  as  1.6  meters  in  this  study,  but  based 
on  other  studies  in  the  area,  can  vary  substantially.  The  compressional  and 
shear  velocities  are  consistent  across  several  different  inversions.  The  Q  values 
are  from  previous  studies,  and  are  not  necessarily  consistent  with  the  attenu¬ 
ation  coefficients  determined  from  the  data  in  this  study. 


A  variety  of  previously  introduced  techniques  which  primarily  depend  on 
surface  wave  dispersion  to  constrain  shallow  shear  structure  have  been  applied 
to  the  problem  of  site  characterization.  These  techniques  have  been  applied  in 
conjunction  with  standard  refraction  analysis.  Since  the  spatial  sampling  and 
frequency  content  of  the  data  analyzed  were  new  and  unique,  a  number  of  sug¬ 
gested  experimental  revisions  have  been  identified. 

The  first  change  to  be  made  in  any  continuation  of  this  work  would  be  ac¬ 
quisition  to  enhance  the  low  frequency  portion  of  this  data  set.  Data  acquisition 
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with  2  Hz  geophones  would  enable  accurate  inversions  to  a  deeper  depth  than 
was  possible  in  this  study.  Several  additional  changes  in  acquisition  would  be 
necessitated  by  the  use  of  lower  frequency  phones.  Table  3  indicates  that  the 
lower  the  frequency,  the  poorer  or  wider  the  temporal  resolution  would  become. 
As  a  result,  the  events  would  need  to  be  spaced  further  apart  in  time,  which 
would  require  greater  source-receiver  spacings  than  used  in  this  study.  This 
would  add  more  complexity  to  the  study  due  to  increased  scattering  and  lateral 
velocity  variations.  The  SWIG  provided  excellent  data  for  this  study,  but  was 
starting  to  degrade  due  to  lack  of  energy  at  the  farthest  offsets.  This  situation 
could  be  improved  at  larger  distances  by  stacking  multiple  shots. 

One  difficulty  with  the  Betsy  data  was  the  existence  of  an  airblast,  which 
contaminated  the  first  higher  mode.  The  airblast  could  perhaps  be  reduced  by 
berming  the  side  of  the  source  closest  to  the  receivers.  Alternatively,  bermed 
explosive  sources  as  used  by  Simila  (1982)  could  be  used,  especially  if  an  at¬ 
tempt  was  to  be  made  to  increase  the  source-receiver  spacing.  In  this  study 
area,  it  is  not  clear  that  improving  the  higher  mode  Rayleigh  data  would  contrib¬ 
ute  much,  since  the  different  frequency  range  precluded  using  the  higher  mode 
in  the  inversion  process.  The  fundamental  mode  dispersion  analysis  would  be 
made  easier  if  the  events  were  more  widely  spaced  in  time  cue  to  a  greater 
source  offset.  Additionally,  if  the  source  offset  was  increased,  a  longer  time  in¬ 
terval  would  need  to  be  recorded.  The  sample  rate  could  be  decreased  in  this 
case,  because  it  is  unlikely  that  any  usable  data  would  be  available  beyond  50 
Hz. 

One  difficulty  with  the  analysis  of  the  data  in  this  study  was  that  there  ap¬ 
peared  to  be  some  lateral  variations  in  velocity  structure.  One  method  for  deal- 
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ing  with  this  would  be  to  shoot  the  refraction  data  from  both  ends  of  the  linear  ar¬ 
ray  and  process  for  dip.  Two  previous  studies  in  the  area,  Bogaards  (1989)  and 
Grant  (1988)  looked  at  this  problem  and  found  no  evidence  of  dip.  Bogaards 
(1989),  however  did  find  some  evidence  of  azimuthal  variation  of  velocity  struc¬ 
ture.  Additionally,  different  studies  have  determined  a  different  thickness  for  the 
shallowest  layer.  It  seems  likely  that  the  shallowest  layer  varies  gradually  in 
thickness  throughout  the  area.  Any  precise  spatial  technique,  such  as  a  refrac¬ 
tion  survey,  would  yield  a  velocity  structure  in  a  small  area  that  might  be  sub¬ 
stantially  different  from  the  average  structure.  A  surface  wave  technique  would 
tend  to  yield  an  average  structure,  and  in  fact  should  provide  a  better  estimate  of 
the  average  as  the  source-receiver  spacing  increases. 

The  final  area  of  concern  would  be  the  determination  of  Q.  Q  appeared  to 
be  poorly  resolved  in  the  inversion  process.  This  poor  resolution  could  be  dealt 
with  in  two  ways.  The  first  way  would  be  to  determine  Q  from  some  other  out¬ 
side  method,  such  as  was  used  by  Grant  (1988).  Alternatively,  some  effort  could 
be  put  into  determining  the  other  factors  influencing  the  attenuation  (such  as 
scattering)  with  frequency  and  correcting  tor  those. 

The  dispersion  analysis  techniques  used  in  this  study  appeared  fairly  ro¬ 
bust.  It  was  easiest  to  determine  dispersion  with  an  array  of  data.  The  combi¬ 
nation  of  MFA  and  PMF  would  make  It  possible  to  determine  group  velocity  with 
single  phone  measurements  if  necessary. 


APPENDIX 

MULTIPLE  FILTER  ANALYSIS  OF  BETSY  DATA 


The  multiple  filter  analysis  (MFA)  is  a  single  source-receiver  pair  analysis 
technique  that  was  applied  to  both  the  Betsy  and  the  SWIG  data.  The  SWIG 
data  results  were  presented  in  Chapter  3,  section  B.  The  Betsy  results  are  pre¬ 
sented  in  this  appendix.  Each  of  the  plots  contains  up  to  four  peak  amplitude 
points  of  the  envelope  or  instantaneous  amplitude  for  very  narrow  band  Gauss¬ 
ian  filters.  The  peak  amplitudes  from  largest  to  smallest  are  represented  by 
squares,  circles,  triangles  and  plus  signs.  The  corresponding  spectral  ampli¬ 
tudes  are  represented  by  the  same  symbols  in  the  top  portion  of  the  plot.  The 
trace  to  the  far  right  is  the  input  seismic  trace  scaled  in  linear  time.  The  second 
trace  from  the  right  is  the  input  trace  scaled  in  linear  velocity  corresponding  to 
the  group  velocity  scale  at  the  far  left. 

The  data  is  presented  in  approximately  5  meter  increments.  The  primary 
features  that  can  be  distinguished  are  the  fundamental  and  first  higher  mode 
Rayleigh  waves.  The  higher  mode  Rayleigh  wave  blends  into  the  airblast  at  300 
m/s  for  the  high  frequencies. 
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Figure  126.  MFT  analysis  at  5.0  meters  for  B6ts>  data 


5  00  4  00 


AMPLITUDE  (CM-SEC) 


AMPLITUDE  (CM-SEC) 


Figure  143.  MFT  analysis  at  80.0  meters  for  Betsy  data 
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Figure  146.  MFT  analysis  at  94.5  meters  for  Betsy  data 
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